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I. PROBLEM OVERVIEW AND DOMAIN DESCRIPTION 


The problem we consider is where an acoustic pressure wave emitted by 
an active sonar impinges on a submarine. In our case this is a double-hulled 
submarine. Active sonar works on the principle of detecting the reflection off 
a solid object of an acoustic pressure wave emitted by a source, by which one 
can calculate distance and direction to that object as in Figure 1, but instead of 
concentrating on the reflected wave we look at the interaction between the 
structure and the incoming wave. This generates scattered pressure waves. 
The scattered pressure waves include waves which decay as they travel 
through the fluid (evanescent modes), and waves which do not (propagating 
modes). Since propagating modes do not decay they can be detected. The 
main thrust of this thesis is to determine the characteristics of the propagating 
modes, such as amplitude and energy. The optimum situation would be to 
perturb the double-hulled structure at a resonance frequency which would 
increase the amplitude of the propagating modes making them easier to 
detect. We investigate the steady state behavior of the propagating modes and 
the resultant shear strain field in the I-beam. At steady state the 
characteristics of the propagating modes stabilize which might be used as an 
acoustic signature of the structure. In addition when the shear strain field 
reaches steady state its value throughout the solid may be used to isolate areas 
of the I-beam where large stresses occur. This information could be useful in 


the design of double-hulled vessels. 


LD) 


Figure 1. Physical Situation 


SUBMARINE 





To reduce the model to one where we can apply numerical techniques to 


simulate its behavior we must make simplifying assumptions which are 


discussed below. 


A. ASSUMPTIONS 


The source of the pressure waves is far enough away so that the 
impinging wave can be approximated by a normally incident plane 
wave. 


There are no other sources present such as might occur with bottom 
bounce, surface reflection and the like. 


The area of interest—AOI as shown in Figure 2 where the wave 
impinges is where the inner and outer hulls are joined or connected by 
a supporting spar forming an I-beam shaped domain. 


The I beam shaped domain is a uniformly continuous linear isotropic 
elastic medium with no cracks, welds or other deformities. 


The dimensions of our AOI are such that any curvature of the surfaces 
can be ignored. 


The center spar or support beam occurs at regular intervals through the 
structure as shown in Figure 3 allowing us to truncate the domain to 
the left and right of center spar using periodic boundary conditions. 


The incident wave does not displace the submarine. 


The cavities A and B in Figure 2 enclosed by the inner and outer hulls, 
and the center spar is void and contain no sources. 


We are only interested in displacement in the x and y directions as 
given in Figure 2 which reduces our problem to two dimensions. 


The fluid is seawater, the solid is steel. 


Our model is now reduced to one where we have two coupled media—fluid 


and solid, and the characteristics of each will now be discussed in turn. 





Figure 2. Area of Investigation 


A. THESOLID 

A cross-section of the double hull of length 8a is shown in Figure 3. Note 
a is a scaling constant for distance. The domain is essentially infinite in 
extent in the x (length) and z (depth) directions. Since it is of uniform shape 
in the z direction, we can neglect this coordinate and reduce our problem to 
two dimensions. Our area of investigation is outlined in Figure 3, and we are 
now faced with the problem of providing boundary conditions in the x 


domain. 
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Figure 3. Cross-section of Submarine Hull 


We have a normally incident plane wave impinging over one boundary 
of the domain. Thus the pressure is the same at all points (there is no phase 


shift) along the fluid/structure boundary. Due to the periodic nature of the 


structure we expect the solid to behave the same over given intervals of 2a, 


that 1s 


f(x,y,t)= f(x+2aN,y,t) (1) 


where N is an integer and f a function describing the state of the solid such as 
displacement, velocity and so on. Using this periodicity we can truncate our 


domain at a and -a and use the following periodic boundary conditions 


f(-a,y,t) = f(a,y,t) (2) 
and 
f*(-a,y,t)= f*(ay,t) (3) 


where k is a positive integer and can signify the derivative with respect to x, y 
or t, depending on the function f. 

The solid has now been reduced to a two-dimensional linear isotropic 
elastic medium whose governing equation of motion for points interior to 
the domain is given by the plane strain elastic wave equation which is 


Ht) a ct uate —— =p - 
ax? dy” ax? oxdy) “° at? 


2 2 2 Zz 2 
1 Soo Zea ou po Xap, oe (5) 
oy x 





(4) 


u and v are displacements in the lateral (x) and transverse (y) direction, 4 and 
A are Lamé constants, and ps is the density of the solid. The I-beam is 
deformed due to the imposed fluid pressure. Balancing normal forces at the 


fluid/solid interface we obtain the boundary condition 


0 
Tyy = a{ 24 4 | +2u aD = poe (6) 


where ptotal is the total pressure, and ty, the normal stress component. All 
other boundary conditions are either traction free or periodic and are 


summarized below. 


Txx = 90 
on surface DH and Elin Figure 4, (7) 
Txy = 0 
Tip = 0 
on CD, EF, GH, IJ, and KLin Figure 4, (8) 
Txy =0 
a total 
on ABin Figure 4, and (9) 
Try = 9 
periodic 
elsewhere. (10) 
conditions 


Txx, Tyy and Tyy are the normal and shear components of stress and are 


represented by 


Ou adv Ov 
= A| —+— |+2u) — J, 11 
wy eae x] ne 
ou Ov Ou 
T xx = {Se a 2u( and Os, 
Txy a {Hs 2) (13) 


Figure 4 is the unscaled domain, in which s, and k are constants used to vary 


the size of the cavities A and B of Figure 2. Note: 0<s<1, 0<k<1. 
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ee ee 
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Figure 4. Unscaled Domain 


B. THE FLUID 
We are interested in the scattered pressure waves generated by 


perturbations at the fluid/solid interface. The partial differential equation 


governing the propagation of these waves in the fluid is the linear two- 


dimensional wave equation given by 


—- = + (14) 


where cf is the speed of sound in the fluid. 

These perturbations arise when a normal plane wave of magnitude P and 
frequency @ impinges on the surface of the solid, causing it to deform. Since 
the pressure wave is of normal incidence (there is no phase shift at points 
along the surface of the solid) and the solid is a periodic structure we can 
expect the behavior of the scattered pressure waves to be the same in given 
intervals of 2a, thus we can use periodic boundary conditions in the same 
manner as was used for the solid. 

Ideally, if the solid were acoustically hard our total pressure would only 
have two components, incident and reflected. In reality the solid is perturbed 
by the incident wave, generating scattered pressure waves which propagate 
out into the fluid domain. The total pressure in the fluid can be represented 


by 


ptotal = pincident ee preflected a pscattered (14) 


where pincident = eT Cag preflected = Bee and pscattered jis to be 
determined. To avoid cavitation at the fluid-solid boundary we employ the 
inviscid form of the Navier-Stokes equation which we refer to as the 


compatibility condition and is given by 


ree ee 2 (15) 


where pf is the density of the fluid. Note that only the scattered term of the 


pressure is included in this equation since 


it and it 


el a = ~ik ge = ik pee (16) 


2p 
oy y=0 


thus 


I R 
ES he | = (). (1 ge. 
oy oy y=0 


The scattered pressure waves generated at the fluid/solid interface are 
composed of propagating (non-decaying) and evanescent (decaying) modes 
which must be allowed to propagate off to infinity in the positive y direction. 
To do this we must employ a non-local radiation boundary condition whose 
implementation will be discussed in greater detail in the discretization 


section of this paper. 
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Il. NON-DIMENSIONALIZATION OF THE GOVERNING EQUATIONS 
OF MOTION AND BOUNDARY CONDITIONS 


The problem addressed models the propagation of a scattered pressure 
wave in two dimensions. This is described by the two-dimensional scalar 


wave equation 


Be 9 SE (18) 


where p is the pressure and cy the speed of sound in the fluid. To facilitate 
implementation and to free ourselves from the requirement of using a given 
system of measurement such as metric or imperial we scale or non- 


dimensionalize Equation 18 as follows. Let 


= — 4+ (19) 


represent the unscaled wave equation. We now use the following 


relationships: 


ie x = — y= p=". (20) 


Here @ is the scaling constant for time and the frequency of the incident plane 
wave, a is the half length of the I-beam and the scaling constant for distance 
and P is the scaling constant for pressure. Note that when taking derivatives 


with respect to x we get 


11 


2 2 
dildo, @ 1d 


as 21 
dx adx dx? a* dx? 7” 


This also holds true for derivatives with respect to ¥, and a similar relation 
results when taking derivatives with respect to t. With the above 


relationships, Equation 19 is now written as 


—. (22) 








Dip Daa 2 Z 
Se (23) 
cf ot’ ax" oy 
Defining 
wa 
kf = Cf j (24) 
Equation 23 is now written as 
2 2 2 
7 a a (25) 
ot” ax” ody 
The elastic wave equation, which is a vector wave equation is given by 
a7u du a°u  a*v a7u 
—~ + — 1+ (A +p) — + = P.—> (26) 
(os =] ( ae aa Psay2 
a*v , dv a*v nu d*v 
oe ee Ae —— = 0.-——=. (27) 
ape =] ( ee =| Ps 542 
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The following relationships allow us to write Equations 26 and 27 in a more 
convenient form, 


P=c2 and ee 


(28 
Ps Ps 


The constants cy and cry are the lateral and transverse velocities of the solid. 


The unscaled forms of Equations 26 and 27 are now written as 








2— 2— 2— 2— 2— 
o|0°u Ocul > O\Nsemne Oo Ou 
1 | ee Pe ae 29 
aS = («t 4) 2 Zz OF? i 
a= 2= a= 2 2= 
DIG 0 mOmD ee \| Gao). Gat Ov 
SA gl = — =-——. 30 
a =D + +(e 4) ae + an ry (30) 


We use the same scaling relationships as before for x, y and t (Equation 20) as 


well as 
a ee” ; 


which are the scaling relationships for the displacement, and rewrite 


Equations 29 and 30 as 








2 2 2 2 2 
2|| (Oo Deiat Dad Oa 7 0°U 
Cr) = + se [ticp — cr |] SMe ts = Dw* —. (33) 
aE: Oe = | : A dy? a* dxdy at? 
Defining 
ee = ie and kr ae (34) 
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and dividing Equations 32 and 33 through by w2 and cancelling common 


factors they reduce to 


1 | of ug geet) (fle 
ke ax? ay" k? kA 


(ato ato) (4 
k2\ ox? oy? } (ke kA 


1 


1 


Chey 
ax? 


Q*v 
—— + 
ay? Oxdy 


Collecting like terms gives us the final form 


tatu 1 eu, 
kf dx? =k dy? 


1, 1%, 
Ke oxt Key 


oe 
kt kt 


oes 
ke kp 


1 


1 


| 


lize 





O-u 
oxoy 





d*v a a7u 
oxdy 


_ ou 
ore 


_ oo 
Ore 


The surfaces in contact with free space are stress free. 


KS) 


(36) 


(37) 


(38) 


(The fluidy selina 


boundary is dealt with separately). Therefore the stresses tx; and tyy on 


surfaces EI and DH and Tyy and Txy ON surfaces CD, EF, GH, IJ and KL of 


Figure > are Zero. 


These components of the stress tensor can be written 


du dv = 
Tyy =H ayox = 
du dv Ou 
ye Se +o ao 
a {Z ar ay 
du dv Ov 
| Se 0 
vw E a) “a 
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(39) 


(40) 


(41) 
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x=-] x= (0) x=1 





Figure 5. Scaled Domain 


where x, y, W@ and 0 are the unscaled components of distance and 
displacement and yp and A are the Lamé constants. We will now scale 
Equations 39 through 41 in turn. 

For Equation 39 using the same scaling relationships as before (see 


Equations 20 and 31) gives 


iS 


aody aox 


Cancelling common factors reduces it to 


cig 


For Equation 40 we first divide through by ps; the density of the solid, and 


using the same scaling relationships we get 


aox ady}) ps a ox 


We know that 
Me 
L/ Ps = Cr. (45) 
Similarly it can be shown that 


= Sr (46) 


Using Equations 45 and 46, substituting into Equation 44, and cancelling 


common factors gives 


y) 9) du dav 2 OU 

=) —+— |+2c7 — = 0 47 
(ct 4)# 4 aay sad 

Or 

Zz 

Ch du dv ou 

——2 | —+— 1+ 2— = 0. “ 
E (x =| Ox = 


Using the previous definitions of ky and kr it can be shown that 
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=. (49) 
ct kt 
and when used in Equation 48 it reduces to our final form 
2 fi 
os az, Ou _o. (50) 
kp ox (kp Jey 
Following the same procedure for Equation 41 its final form is 
2 2 
Hf _|u , Hav is 
kp Ox kf oy 


At the interface the fluid cannot exert a force tangential to the boundary, 


hence the shear component of stress is zero there, and is given by 
ou ov 
Try = BM) += | =0.- 2 
—* e a ee) 


For the normal component of stress at the interface we refer back to Equation 
6, Section B and Equation 15 in Section C of Chapter I to see that the normal 
component of stress is given by 
ou ov ou Res 

Ty, = Al —+— |+2u—=-(p +p +p? |, 53 

; e =| ps =p! + p® + p*) (53) 
where the superscripts I, R and S signify incident, reflected and scattered 
pressures. The incident and reflected pressures are given by 


1 _ -ikpy-it ik¢y-it 


p and ps =e 


Txy is scaled the same as before. To scale tyy properly we will need the 


compatibility condition, whose unscaled form is 


i 


Pf ae2 (54) 





We refer back to Equation 17, Chapter I, Section C to see that only the 
scattered pressure need be considered here. 


Scaling Equation 54 we get 





P op” Pla 
a = ppo*D—, (55) 
Or 
2 
p° Diet 
—~P = pw aD—~. (56) 
oy y= Q “i at? 





Equation 56 gives us a convenient choice for the scaling constant for pressure 


of 


a 
P=@ aDpy (57) 


and when substituted back into Equation 56 cancels as a common factor to 


give 


=—. (58) 





Dividing Equation 53 by ps and scaling gives 


(59) 





AD{ du dv) 2NDdv_-P 
20 (2 = }: Ll —(p! + p® +p’) 


p.alax dy) pea oy ps 





ox dy UF 


Substituting the value of P from Equation 57 we obtain 
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AP (Su, 22), HP. 2 = PL y Papi +p +p°) (60) 
psa\ox oy) psa OY ps = 
or 
ky Ou dv ov Die eR OS 
— — + — 2_—--— k + 61 
Ei E =| 0 ekt(p ae = ia 





where € = pf/ ps. 


Evaluating the incident and reflected pressures at y = 0 we get 


pl = pR = eit, and when substituted into Equation 61 gives 


ke du dv ov 2 Sie es 
— —- 24 —+— |+2— = -2ek — Ekfp”. 
E | ae =| ay Se ine (62) 


uy, 


II. DERIVATION OF THE FINITE DIFFERENCE FORMULAE FOR THE 
GOVERNING EQUATIONS 


Throughout the derivation of the finite difference approximations we 
identify gridpoints of the scaled domain with the subscripts i andj and 
superscript n. Variation in the x direction is denoted by 1,1 <i<N, where N 
is the total number of subdivisions (since we are using a square grid the 
number of subdivisions in the y direction is the same), variation in y with j, 1 
<j<WN, and time with 1,0<n<eo. Lower case indices denote varying 
quantities, while upper case letters denote fixed quantities. ere would be the 
value of the scattered pressure for all values of 7 at the grid level y = KAy, thus 
Dh is a vector, while Pi is an element. 

As was discussed above we use the letter 7 to denote variation in the x 
direction. This is a common convention and we do not wish to deviate from 
it. To avoid confusion with the complex quantity, ~-1, also commonly 
denoted by i, we state the following rule, that whenever the letter 1 appears as 
a superscript it denotes the complex quantity /-1, and when appearing as a 


subscript it is an index denoting variation in the x direction. 


A. FINITE DIFFERENCE APPROXIMATIONS FOR THE EQUATIONS 
GOVERNING THE BEHAVIOR OF THE FLUID 
The scaled domain of consideration as shown in Figure 6 has periodic 
boundary conditions applied at x = 1 and -1 and a radiation boundary 


condition for the propagating modes at y = 2. To model the two-dimensional 


20 





wave equation numerically we must derive an equivalent finite difference 
formula, from which we can solve for the pressure for all subsequent time 


levels as follows: 


ke = + (63) 


is the two-dimensional wave equation as derived in Chapter II Section A. 
Using central difference approximation for all second order partial 


derivatives the equivalent finite difference equation for Equation 63 is 


isa n+1 1 n n n i n n n 
Palla ia 2p; + Pig = 5 [Phen j ~ 2Phj + Pha j|+ col Ph jaa ~ 20 + Phy 


(64) 


where h = Ax = Ay is the step size and Aft is the increment in time. The 


truncation error for Equation 64 is O(h2) in space and O(Af2) in time. Solving 


for py explicitly we have 


2 2 
ip -[2- 2 ss 2 51 Pha,j + pe 113) ie eel eee 1 =p. (65) 
kfh keh 


Z 
Letting p* = a Equation 65 can be written as 


n+1 _ 


Zeer a sass n n n n—1 
oan 2(1-2p \pr fe pj Be iS * Pi i+ al Sigs (66) 
The Von Neumann stability criterion 


(67) 


aI 


C= 


must be satisfied to ensure the stability of Equation 66.! 

Special attention must be paid when applying the wave equation along 
the boundaries at x = 1 and —1 for it is here that we make use of the periodic 
boundary conditions. Applying Equation 66 at x =-1,i=0 and y = jh (see 
Figure 7) yields 


n+1 _ 


2 2 —j 
po; = 2{1-2p" )pp,;-P ee Pj + PO, jot + PO, — Po, j (68) 


This requires the value of p at the point (-1, jh, nAt) which les outside the 


domain. By using the periodic boundary conditions 


pt, y, t) = p(-1, y, #) (69) 
op _ oP 
ox |(+Lyt) ax |(-Lyt)’ (70) 








we can substitute the value of ((N-1)h, jh, nAt) (where N is the total number 
of subdivisions in the x direction) for the value at (-1,jh, nAt), allowing us to 


evaluate the wave equation at the boundary. 


B. APPLICATION OF THE RADIATION BOUNDARY CONDITION IN THE 
NUMERICAL SCHEME 


As was mentioned at the end of Chapter I (Section C), we apply a non- 
local radiation boundary condition (referred to as nlrb) to the fluid to 
simulate an infinite domain in the positive y direction. Our domain is 
truncated at x=1 and -1, forcing our fluid domain to act as a waveguide. The 
scattered pressure can be represented as a series of plane waves which take the 


form 


1For treatment of the von Neumann Stability Criterion see Appendix A. 


Zo 





p ia 


2 
ey 
O 


kj - Vp for propagating modes 


elrextBy-!) where Be (71) 


il vp - kj for evanescent modes 
and 7 = kz. We note that when fj = in] Vz ~ kf Equation 71 yields 


o Brytil Yex-t) 


which is an exponentially decaying quantity for positive values of y. This fact 
allows us to neglect evanescent modes when applying the radiation boundary 
condition, at y = 2. 
The total scattered pressure is given by 
penne aye EEA) (72) 
k=—00 
This is composed of propagating and evanescent modes. Far from the fluid 
solid interface where only the propagating modes are assumed present (for 
reasons given above) the scattered pressure is written 
M 
p — SS ayet ee Bay-#) (73) 
k=-M 
where M is the total number of modes (positive or negative) under 


consideration. At the boundary y = 2, we apply the nlrb operator (Scandrett 


and Kriegsmann, 1992, unpublished paper). 


Op op 
B,(p)=—- + By (74) 


to the individual modes of the scattered pressure of Equation 73 which yields 


J, 


M 


B(p)= > sane Veet Pky! 4 »> bil 5 -{ ape? Ykx+ PRY )) (75) 


k=-M 


Equation 75 reduces to 


M 
B(p)=  d_(iBy — iB Jay ef 142+ Bayt) (76) 
k=-M 


The right-hand side of Equation 76 is identically zero. The boundary operator 
has annihilated the propagating modes and since the evanescent modes are 
assumed to have negligible magnitude there, any scattered pressure waves 
reaching the boundary experience no reflection, simulating an infinite 
domain in the positive y direction. 


To apply the nlrb operator at the boundary y = 2, we rewrite Equation 75 as 


-> By <(a eT) Belh+ 742) ) = 0 (77) 





where we evaluate the pressure at a constant time T and along the boundary 


y = Jh (i.e. at y = 2). We incorporate ag and e~7 as aye“? and define this to be a 


=1 for all valuecuara 


new constant oj(T). Note that |a,(T)|=l}] since fe“ | 


ax(T) is unknown so we must derive an alternative expression to be able to 


evaluate Equation 77. The k‘h propagating mode can be written as 


Dy = ayetl Yk*+Bry-t) (78) 


and when evaluated at y=Jh and at t = T and employing our new constant 


ou(T) we obtain 
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(T)e@PKIh gt kX, (79) 


To isolate a, (T)e'Palh we multiply Equation 79 by eYkX (we apply 


orthogonality) and integrate over the x domain to obtain 


1 
| p(x, Jh, Tet’ dx = 20% (T)e!Pal? (80) 
~—] 


OT 


1 
or, (T)eiBaIh = ; { p(E,Jh, Tet dé, (81) 
=| 
where ¢ is a dummy variable of integration. We substitute this value of 
a. (Te? KIF back into Equation 77 to obtain 
Op 1 Be iy, (x- 
he + jh, Tek -S)de = 0 82 
fil am 
which allows us to apply the nirb at the boundary y = 2 and from Equation 81 
we will be able to evaluate the amplitudes of the propagating modes of the 
scattered pressure. 


We now proceed with deriving the finite difference approximation for 


the radiation boundary condition. Using a central difference approximation 


for ue _ 7,7 Equation 82 is written as 
ay 


27. 


Poe eee 1 SB, 2 (E, Jh, T)ebkC-S)de = 0. (83) 
2h 2 oy ae 


The trapezoidal rule for integration is 


ff(x)dx = =(f + 2 fo + 2 ie Bs | + i) (84) 


a 


and is used for the integral in Equation 83, however it can be more compactly 


expressed as 


r=) one 


Nie 


i 
[flxWdx=hD5,f, where 6,= (85) 


r= 


— 


elsewhere, 


where / is the total number of nodes in the x direction. When substituted 


O 
into Equation 83, using a central difference approximation for (C fia e) 


Equation 83 can be written 


Pipa ~ Pi j-1 ee i = 3 1 
una Bele §,efte(ti-$r)( yn 0. 86 
2h ty Pe Hy ; aie ee 


Lae 2At 
Multiplying by ]—, we have 


2A . 
a (p? je +d Das kOr ca KI (p mat prt) =0. oe 
r= k= 


Define A to be a matrix whose 1,r entry is given by 


M . 
r)= S By.6,e1% 4~ Sr), (88) 
k=-M 


Upon substitution into Equation 87 we obtain 
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2At 
7 (Phys Piy- 1) + Apps — App y* = 7 (89) 


n n al n= 

We note that PiJ+i’ PiJ-1 Pr J and Pry are all vectors due to the J] 
subscript as described at the beginning of this chapter. The radiation 
boundary condition and the wave equation must both be satisfied at the 
artificial boundary y = 2. Applying either of the two conditions will require 
the use of pseudonodes which lie outside the domain. Through the 
combination of the two equations we will be able to eliminate this 
requirement. Reproducing the two-dimensional wave equation and the 
radiation boundary condition, (substitute k for all x indices in the wave 
equation and radiation boundary condition, since these are dummy indices), 


we have 
n+1 g) n n—1 At? n n n n 4 n — (0 
oan —Pk,} * Pk] Re Peet Phe, Phys t Py a ae (90) 


2At 1 
12 <> (Pe pain Pk. J- 1} + Apes} — Ape = 0. (91) 


Note that both equations are being evaluated at y = 2 and that the vectors 
Prt, J and Py +1] have a circular shift and are of the same dimension as the 


rest of the vector in Equation 91. Collecting like terms in Equations 90 and 91 


we obtain 


‘ “50 Pk tof + Pky + PE ge 5 |Phay + Phas, ~ 47h |- 20k) = 9 
KR? k,J+1 — kPh? Syl ey a kn? k-1,] © ?k+1,J k,] k,] 


(92) 
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2At 2At 
Ne i oe Ae oe 1 eaeen hig ee ae = 0. (93) 


At 
Adding —7; times Equation 93 to Equation 92 yields 
2k 


f 
At? re At 1 At 
arty ore eee 


4At At? At? 
+| ~~ - 2 |p? ea eae) 
| KPI he ik Kh 5 Pk Yo KP 12 Seen (94) 


Solving for Pe, in Equation 94 we obtain 
7% 2 2 
[+—yA +| 2-——> 


Ae zs poe zy oe 


The terms differenced in x, that is Pky / By Jv Dre y can be more compactly 
4Ate 
272 
At k¢h 
242 on the sub and super diagonals. We must also allow for 


expressed as TP, where T is a tridiagonal matrix with 2- on the main 





diagonal and 





the periodic boundary conditions when constructing T. To do this we replace 
Z 

the (N, 2) and the (1, N-1) elements of T with . 5 

matrix. The general form of T can be seen from Equation 16 Appendix C. 





where T is an N by N 


Equation 95 is now written as 
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-1 
2 
n+] At 2At” py n At i=l 
Sioa | | Sa e +P 
Pk,] | 2 | Ph PkJ-1* “Pk, ] 2kP Pk,] (96) 
which satisfies both the wave equation and the radiation condition at the 


boundary. 


C FINITE DIFFERENCE APPROXIMATION FOR THE ELASTIC WAVE 

EQUATION 

To investigate the propagation of disturbances in an “I-beam”-shaped 
domain as shown in Figure 5, it will be necessary to apply the elastic wave 
equation to points in the interior of the domain. (The boundaries will be 
dealt with in a separate section.) By only considering displacements in the x 
(lateral) and y (transverse) directions, the problem becomes one of plane 
strain in two dimensions. Reproducing the scaled equations derived earlier 


for motion in the lateral and transverse directions 





1d¢u 1 atu 1 1)\a4 dtu 
ae ae eT em ee) 2 
Ki Ox kr oy ky kr Oxdy ot 
2 Z 2 Z 
wn een ES 8 
kr ox ky oy ky kt oxoy ot 


we use central difference formulas for the partial derivatives in Equation 97 


to get the equivalent finite difference equation 


1 n n n 1 n n n 
nk? eens 4 2u; ote uty == atin = 2u; = el 
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1 ae n n n 
oF . +] Wy (oF jaa Ci-1,j41— P841,;-1 7 44-1, j-1 


9 n+1 is! 
= —|ult) — Qu ult (99) 


The truncation error for Equation 99 is O(h2) in space and O(At2) in time. 


: 1] 
Solving for Wij 


explicitly in Equation 99 we obtain 
+2 


2 
ae At ut n A ut n 
u. Pe ea as oe 
keh? 1,j+1 eal 


tu: 4. 
oe eR eeciey i 1 


2 
At i 1 | n n n n 
Sa a a a Ue . pn of ae . — WU. Sey tad ae 
Ah? s "7 +1,j;+1 al pas Ll a bf 


ZAte | ot il _] 
fe 73 aan 2 |i (100) 


Using the same method for Equation 98 we solve for v; * to obtain 


; 
ott = (ot +ur J+ AE (on +r 
1,] ken? rly ine keh? 7] nt | 


2 
+—} —- -—- —— (u! ee eee ee ae 1} 
mee a 21+ ie — 1) ee ee — lf ol 


Zaye 1 = 
2-3 h? sta ae ye 


To ensure the stability of Equations 100 and 101 the Von Neumann stability 


condition of 


oe 


Ae <= (102) 


must be satisfied.? 


D. APPLICATION OF THE STRESS-FREE BOUNDARY CONDITIONS TO 
THE SOLID 
The boundary conditions of the two-dimensional domain are broken 


down into two major categories, those whose normal vector is (( , where 
Txx = Txy = 0, and those whose normal vector is a where Tyy = Txy = 0. 


cto 
These are in turn divided into two classes. For n = ( Jthey are 


1 
al. The boundary whose unit normal is a | , that is facing in the positive x- 


direction, the surface EI in Figure 5 and 


-1 
a2. The boundary whose unit normal is (( , facing in the negative x 


direction, the surface DH in Figure 5. 
0 
Similarly for 7 = 1} 
0 
bl. Boundary whose unit normal is | , the surfaces AB, GH and JJ in 
Figure 5 and 


0 
b2. Boundary whose unit normal is A.) , the surfaces CD, EF and KL in 
Figure 5. 


The application of the stress-free boundary conditions for cases al and b1 is 


discussed below. 


2For a brief treatment of the Von Neumann stability criteria see Appendix 
B. 
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1. Application of the Stress-free Boundary Condition for the Case of 


i 


The boundary under investigation is identified in Figure 8 as XY. 


(N-1,/+1) 





1 
Figure 8. Boundary with Normal 4 | 


The governing equations as derived in Chapter II Section A 


du Ov _ 


AME SoG 103 
ey) ay Ox a 
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Z Z 
on-(G)et (ih a} a 





k? ox ky y 
2 2 Z Z 
cS ae aes av aoe (105) 
en Ox kr oy ky kv oxdy ot 
2 2 Z 2 
oR A aos pees (106) 
kr ox ky oy ky kr Oxoy ot 


Concentrating on displacements in the lateral (x) direction, a typical boundary 
point as shown in Figure 8 must satisfy the boundary conditions, Equations 
103 and 104 and the governing equations of motion, Equations 105 and 106. 

If we apply Equations 105 and 106 at the node (N,;) in Figure 8 we will 
require values for Wy 4q i/0NG1 j41/ON41,j-17 N41, j/UNGL je AMF UN jv 
which lie outside of the domain and are called pseudonodes. To eliminate 
this reliance the technique as developed by Ilan and Lowenthal (Ilan, 1976, pp. 
431-453) is followed, and is presented here. 

The lateral displacement (u) at the node (N-1,/) in Figure 8 can be 


expanded in a Taylor series as 


higher 

DY eG 5 

ut sul - =) a “4 | meren. (107) 
J Ox N,jn 2 | ox N,jn terms 


where UN-1,j denotes the value of u at the node (N-1,/) and at time level n. 


zu 
Alternate expressions for 5 and 5x2 are given by equations 104 and 105, we 


have from Equation 104 


oo 


ou ke (ke 0 lav (ke av 
a; are ye Pin 2-5-1 awe (108) 
kp \ ky D kt y 





9° a°u kf d° © one. 
wie Swine rei! ia ) 
Ox dt” kf oy kp ke jdxdy 
—9du kp atu (ke _)0*0 
=i ee ape . (110) 
at” ke ay* \ke — jaxdy 


Thus the lower order terms of Equation 107 can be written as 


2 2 Pa 2 2 a 
uN ag = YN aa duane a a “Ly 7?) aa 
kr Cy ot” okpoy’ lk oxdy 





dv od42u d2u 
Using centered differences for Oy Pe): ay2 and the following difference 


formula for the cross derivative term 


07y il 


pkey ai Ti 
dxoy 2h i 


n n n n 
(of, j41 (ON Ne ca 
the finite difference approximation for Equation 111 is 


n ee h 2k? 1 1 n n h* kf n+1 n n—-1 
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kr 
h7ke 1¢ 2 ee 
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Cancelling common factors Equation 113 reduces to 


n eee 7 1 kf n n h?k? n+] n n—1 
Oe; ~ "Nj *| 5742 (ON jaa PN, ja} +55 5 UY = 2a, +N] 


ey. 
2 ke 
ae eel 
=a 2 + || Ee - = fon ON pa OR 4 1} 
sr Niet uN j UN j-1 # a \ON,j+1 ON, j-1~ °N-1,j41 7 N-1,j-1 
(114) 
Solving for explicitly we obtain 
2 
eal 2 Ata) wal “4 n-1, ote 
WG | a eee | KEN 
Ae aye ON +ur anes us (OR, —vU\ 
2 
ae ee a Gils) 


The truncation error for Equation 115 is O(h3) (Man, 1976, pp. 431-453). Using 
the same procedure for the transverse displacement an explicit expression 


n+1 . ; 
6) N,j 1S given as 


Ag Ai | ae 1 2 
<5 (on, m1 ON Ge 1}+ ke ie (uy =| = =the ia) 


Wk 2h 
At ——— (ut agi | (116) 
372 k? Ka N,j+1 N,j-1]' 


a7 


2. Application of the Stress Free Boundary Condition for the Case of 


iC 


In the case of b2, the governing equations 


lov 
ke ax? 


Ou dv 


—+——=(0 
ye oy ox 


; _(kp_, lou, kf dv _ 
Yk ex Tr oy 


1 du 


id'u, Tou fd Lae cam 
kt ox’ ke oy? | kp z adxdy ar? 


Lav f 1 ie ua 
op a | oS 
ke ay* \kP “2 axdy dare 


(117) 


(118) 


(119) 


(120) 


and a typical boundary point is depicted in Figure 9. Expanding in the vertical 


direction at the node (i, M+1) and at time level n, and ignoring higher order 


terms 


it = 4 +i | ali au 
1,M+1 1,M oy — ay? 


using the substitution 


from Equation 117 and 


O*u 
ay? 


Z Dae 2 Ze 
ere pei d e 
k oxdy 








di kee 
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1,M,n 


(121) 


(122) 


(123) 


from Equation 119, Equation 120 now becomes 





n2k2 a2u hh? kk O2u hh? ( k? Q02v 
er ee eee | 1 124 
Ui M+1~ 41M ( =) 2 ao 2 k? ax? 2 kf oxdy as 


(i+1,M) 


FREE SPACE 





0 
Figure 9. Boundary with Normal Oy) 


| diff i er dee th 
We use central difference formulas for jo ae and for the cross 


Q2 
derivative term dy the following finite difference formula is used 


du Ps 1 nN Eaaatt ait n (125) 
dxdy 2h2 Gary ESL “ATi “Ayes 


The analogous finite difference equation is now 


Diep 
1 h“k = 

n 23) n ante fl it) eon n-1 

Ui M+1— 4m is cae ofa) too ee (urN 2u; yt ua) 
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ye eo 
: (ur — 2ui y+ a) 


2 eae 
h? ( k# oer , n n 
A(t : |e (OPM ~ P41,M ~ %-1,M41 7 vPaM+1) (126) 


Cancelling common factors yields 


1 

n n n n 

u; =Uu- -=(v} ee Jaa u 
1,M+1 1,M D) i+1,M i-1,M 9 142 


fu —2u?,,+ur \+ 1_ kp (or +Ur 4 4,-U: +Ur 
2k? i+1,M 1,M i-1,M 4 Ak? i-1,M 1+1,M i-1,M+1 i+1,M+1 
(127) 
1 oo 
and solving for Beh explicitly, we obtain 
2 2 
nad 2 Ate eal 7 i=). Saree 
Us, =| 2—-—s—| St [ay Uy ge He; 
1,M | hh? 3 ké M,j M,] h*ké 1,M+1 
2 Z 
a y2 ag (fea i Mam) 5.3 3 - + |e = ts 
Me jos (v" _ yf } (128) 
oh k? ka i+1,M i-1,M 


eae se Ik 
Similarly for v;,, we have 


2 2 
UVM — 2-26 UM %MT hk KP 0; M41 


h? \ kp ke 
2, 2 
Ae os pi ie | mieel ie ; 
+—-—>~|V; + Ve — | oe oT id) —U; 
OE | i+1,M l 1M} (a "| i-1,M+1 ulm) 
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E. FINITE DIFFERENCE APPROXIMATIONS FOR THE CORNER NODES 

The I beam shaped domain has four 270° corners which are identified in 
Figure 10 as 1, 2, 3 and 4. The treatment of these corners falls into two 
categories, a) corners 1,3, and b) corners 2,4. Within each category the finite 
difference formula applied is identical, the difference between them comes in 
the cross derivative terms of the elastic wave equation. Each category will be 
discussed below. It is important to note that only the governing equations of 
motion are applied. The stress free boundary conditions are omitted due to 
the complexity that arises in trying to apply stress conditions at the corner 
node. We assume as in Fuyuki and Matsumoto that the consequences of 
neglecting these boundary conditions is minimal. (Fuyuki, 1980, pp. 2051- 
2069) 

Category a (Corners 1,3): An arbitrary numerical mesh with Ax = Ay =h 


about corner 1 is presented in Figure 11. The governing equations of motion 





are 
10%u 10% (1 1)\0a% 2 
ce ee ee (130) 
ky Ox kr oy ky kr oxoy ot 
2 2 2 
ak cee oe pis ou = a (131) 
kr Ox kj oy ky kr oxoy ot 
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Figure 10. The Corner Nodes 
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At time level 
[=)). 


(,j+1) 


(+1,+1 


(i+1,j-1 





Figure 11. Corner Node—Category a 


In Equation 124 which describes lateral motion, the normal central difference 
d2u d2u d4u 


formulas are used for the a ae ay?’ terms and for the cross derivative 
d2v 
term Oxy the difference formula of Fuyuki and Matsumoto (Fuyuki, 1980, 
pp. 2051-2069) is applied at the node (1,j) which is 
oa 


53a odie D*D? + D*D" |. (132) 





Where Dy is the forward difference formula 
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p* = uta, j-u | (133) 


and D* is the backward difference formula 


1 
a. n 
DE = [uj wea jh a 


The resultant difference formula for the cross derivative term in v at node i,j 


is 
0*u Nie 7 F P ( 
oxdy i on? ak Yi41,j 7 my 2, pa oe le CEL ei ~ ea - 20% (13) 


which is O(h2). The finite difference equivalent of Equation 124 is 


(| nN il n 
al Care 2Uu nul [+ ae hie 2u," a A 


= 8 ieee (o! +0 no +ue —v. - 2vf',] 
oh k? ka. oie i,j+1 ° “i,j-1 Oi j4 re bs || Di 


mz ] n+1 n—-1 
cat Si a } (136) 


Solving explicitly for wi, from Equation 130 we obtain 


2 2 Z 
n+] At ( n nN a 2At 5 1 1 =} n At n n | 
a — yee yi SS + — } ly. . + ———| 1. . 7a 
i-1, ? 2 DB 2,24°1,)4+1 1,j/-1 
J k J keh J J 
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ey) kph? i+1,) 1? kf 
Wi 
At il 1 ot nN n n n )- n-1 
coo ——_—emee P oe yn art . . — 7. ° —2Y. U; 
ale k2 ler Cet CH iia Or Vl Pi41,j+1 ee! i,j 1, 
(13%) 


ie 
similarly for Equation 125 an explicit expression for 0; j 3s 


2 2 
v; +0: emma C ot Sse | Ue uy + Oe e 
ay “35 Se) joa 12 kp Ne 1,] keh? oe ee — | 


ate ae Salen +u; Ue a — ur —2u!" tee 
oh k? K2 al 1, ely Me pel i+1,j+1 iit ey ao ‘ 
(138) 


For Category b), the governing equations of motion are the same. The finite 
difference approximation of all terms with the exception of the cross 
derivative are done in the same manner. For the cross derivative term in the 


case of Equation 124, the difference formula applied at the node (1,/) is 


Q° 1 
axdy a i) Dy D? + D*D? (139) 





020 
Thus Oxoy at that node in Figure 12 is represented by 


= 2 |20F + of _ tot, -vl, vt (140) 


n n 
Polgei ” CAS eo) > CRS) ae ea - vf, | 
The finite difference approximation to Equation 124 is now 


il 3 1 
ae py kane \+ ( aan aia, 
h?k¢ ae Mi “1, ne Mi fed 8G, j TG jd 


] ae e n n n n n 


=—5(u n+l -2u" ur" (141) 
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Figure 12. Corner Node—Cateogory b 


1 
Solving explicitly for ui we obtain 


2 2 Z 
n+1 -_ At n n 2At 1 1 n At n n 


At? it 1 Qyft n n n n n n n—-] 
"onl Re Re Us jt Cis j-1 7 F-17417 Yisl,7  2i-1,j 2, j+1 97, j= rr 
Ee (142) 
are . n+1, 
in the same manner an explicit expression for v; j is 
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2 2 
n+1 _ At? S52 Nn 2At ] ] N At nN nN | 
v; -+U- 4 -|+| 2—-—3—-| tas Us tos lO FF 
1,J nek es ly h2 k? ke ay) h2k? jer LOS esi 


At? 


] 1 n n n i i oe 
5 ; 4 | ip t Mien jar t Majo — Yi, oj ja W-1] a 


ij 
(143) 


F. BOUNDARY CONDITIONS AT THE FLUID/SOLID INTERFACE 

The boundary conditions at the fluid/solid interface are modified by the 
introduction of a normal plane wave incident on the surface of the solid. As 
a result the normal component of stress is no longer zero. We must allow for 
the effect of the scattered pressure at the interface, which is a result of the 
compliance of the I-beam structure and reflections of displacement waves 
from internal boundaries. This is done by use of the compatibility condition 


as derived earlier. Our necessary equations are 


= — + —_ =() 144 
4 Ox oy eke 
ké ou ke dv ee 
T..,2|/—~—-2 — + —— = -2ékz ~ ek? 145 
yy if FE k? oy 7 iP ic 
2 2 2 Z 
Bie a (146) 
ky Ox kr oy k? ke oxoy ot 
Z Z ja 2 
pe cle nt a7 
kr Ox ky oy ky k oxdy ot 
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a (148) 
oy f or 
2p. dtp, ip i 
Por x2 ay 


where € = pf/ps and p’ is the scattered pressure along the boundary y = Kh (at 
the interface, see Figure 13), and at time level n. 
n+1 


Explicit expression for He andv:, are derived in the same manner as 


before and they are 
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Figure 13. Fluid-Solid Interface 


Looking at the compatibility condition, Equation 142 we use central difference 


approximations to obtain 


n n 
Pi,K+1~ Pi,kK-1___ 1 ( 


2h Bz EK ~ 22 K+ ox) ye 


assuming that it is applied at the boundary y = Kh in Figure 13. 
Solving for p;',_, we obtain 
2h +1 n-1 
pi g-3 = (tk — 207%, +0; 5 \+ OP Ka (153) 
all quantities on the right hand side are known. We now have a method of 
calculating p;,_, which is required when applying Equation 149 at the 


boundary y = Kh. 


Solving for the p”+! term of Equation 149 we have 
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i 4At At? 
Pik -[2- so aie 1K + Pist,K + Pi,K+1 * Pi,K- i}- pix. (15.4) 


The components of the right-hand side with the exception of p/,_, lie in the 
interior or on the boundary of the fluid domain and can be evaluated, 
however p;,_, is essentially a pseudonode for the fluid and is determined 


from Equation 153 above. By this method we have now generated the 
scattered pressure waves caused by the vibration of the solid at the fluid/solid 


surface. 


G. PROGRAMMING CONSIDERATIONS 
The program for this thesis was written entirely in Matlab 4.0 Beta 


version for two reasons, 
e¢ Ease of programming 


¢ The ability to generate quality graphics. 

Although originally intended as a linear algebra toolkit the above features 
have caused it to be used more and more as a high level programming 
language. In our case Matlab was convenient since the fluid and solid 
domains are square matrices, which are easily manipulated in Matlab. 
Updating values is done somewhat differently than in FORTRAN or C and is 
discussed below. 

Equation 66 of this section updates the interior points of the fluid domain 
and is given by 


n+1 _ 


2\ on 2}. n n n n-1 
i 2(1- 2p \pr, +p Peat Pag + Peja tp aay 4 


an equivalent FORTRAN statement might look like 
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DO 10 I = 2,K 


mO.2Z0 J — 2,K 


EeeieeNael) = 20) ,.0-2,0=(RHO**2))*P(1,J,N) 
& mero ee je (Pe (i—l,J,N)tP(I+1,J3,N) +P (1,59-1,N) +P (1, J+1,N) ) 
& moe, ONL) ) 


20 CONTINUE 


10 CONTINUE 


Here each value is updated individually by using a double DO loop. In 
Matlab this is done by “shifting” a grid the size of the interior around the 
appropriate matrix and weighting terms. The equivalent code in MATLAB 


would be 


PNEW(2:K,2:K) = 2*(1-2*RHO%2)*PCURR(2:K,2:K) 
+(RHO%*2)*(PCURR(1:K-1,2:K)+PCURR(3:K+1,2:K)+PCURR(2:K,1:K-1) 
+PCURR(2:K,3:K+1))—-POLD(2:K,2:K). 


where PNEW contains the new values, PCURR the current values and so on. 
All updating is done in this manner eliminating the requirement for 


multiple do loops. 


oy 


IV. NUMERICAL RESULTS 


In an effort to verify our code we checked the behavior of the fluid and 
employed energy conservation methods to check for the consistency of the 
coupled domain. For the fluid waveguide we want to ensure that the 
propagating modes behave as expected, that is, they should not reflect from 


the artificial boundary. This was done by placing a driving force of the form 


n= A, el(Bnyt Yn2-t) (155) 
ee) Qa 

where £,, = ak? —Yn and ¥,=nz and kf = Cf’ (Ref. Equation 24, Chapter 1, 
Section A) at the boundary y = 0, which excited the fundamental, first and 
second modes (n = 0, 1, 2). The coefficients A, had value 1 for all n. The 
amplitude of each mode was measured at the boundary. As can be seen from 
Figures 14 a, b, and c the fundamental and first mode approach the value 1 

with the second mode showing the same behavior but at a much slower rate. 
To check that the coupling of the two domains was working correctly we 
eliminated the cavities of the I-beam which reduced our problem to one 
which could be solved analytically. For a normally incident plane wave only 
the fundamental mode was excited. Since the incident wave displaces the 
solid only in the y direction v has no x dependence and uw is identically zero. 


The values of Ag and v are given by (Scandrett, 1992, interview), 


and v=e# (cyeltLY + coe KLY | 
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For kf = 1 and ky = 0.2542, c; and c2 were given by 0.0625 + 0.037 and 
0.0585 - 0.0374i respectively. In this case the magnitude of ||Ao|| of 0.1211 
compared favorably with the numerical solution of 0.1255. A major 
discrepancy lay in the numerical and analytical values of v (transverse 
displacement). The average absolute value for the numerical solution was 
13.12 while the average absolute value for the analytical solution was 0.0965. 
They exhibited similar behavior but the numerical solution was translated by 
a constant term. We believe this to be a result of there being no displacement 
term to compensate for the effect of imposing a periodic pressure 
instantaneously in time which has caused our solid domain to drift or 
displace, violating one of our initial assumptions (see Chapter I, Section A). 
To nullify this effect we take the time derivative of our steady state solution 
for v and then compare with our analytic value as can be seen in Figure 15. 
Again the numerical and analytical solutions give close agreement. A second 
check was to apply energy conservation methods to our steady state solution, 
from which it can be shown (Scandrett, 1992) that the propagating modes 


must satisfy 


dx. (156) 





Da 


M 1 ] 5 
> aA? =-28oRe(4o)=—2im fry «oH 
Lal 2 Fy 


The quantities 61a and -2ByRe({Ag) for the various values of ky 


are listed in Table 1. Considerable discrepancies exist throughout, which may 
indicate either an error in the code or the inability of our simulation to model 
the high frequency components which exist at the interface. Another factor 
which may contribute to the failure of the integral is the translation of v 
which was discussed previously. With this in mind we must possibly 


consider the following results as being inaccurate. 
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The first aspect we look at is the amplitude of the propagating modes of 
the scattered pressure. For ky = 7, (this we call our original system) the time 
history of the amplitudes for the fundamental, first and second modes are 
shown in Figures 16, 17, and 18 respectively. These indicate roughly the rate 


of convergence of the numerical method. 


TABLE 1. ENERGY CONSIDERATIONS 


M 2 
din--M Bln 


pos | 4azoi6 | 5205 | 02738 
poo | 4a | 0742 | 97 


To gain insight into the characteristics of the solid, we treated the flange at 



















the fluid-solid interface as a thin plate and tried different values of kg 
corresponding to the resonant frequencies of a thin plate with prescribed 
boundary conditions. The first case was a plate with the fixed (referred to as 


FF) boundary conditions given by 


v(—-1) = vO vy (= 
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The resonant frequencies are 


crhn?n* 
Gai) > 
(Lalanne, 1982, p. 103) where v is Poisson’s ratio, h is the thickness of the 
plate, cr is the transverse velocity of the solid and a is a scaling constant for 
distance. 
For the second case, clamped (referred to as CC) boundary conditions were 
prescribed, which are given by 


Veo vl) = 7 el) =) = 0. 
The resonant frequencies of this system are given by 


wo, =p (158) 
4a,/6(1- v) 
where X? = 22.37, X4 =61.67 and X$ =120.9. (Lalanne, 1982, p. 103) 

The purpose of this experiment is to determine whether the I-beam 
exhibits behavior similar to a clamped or fixed plate, since it would be 
advantageous (numerically) if we could substitute a periodically placed 
boundary condition for the remaining portion (lower flange and center spar) 
of the I-beam rather than having to calculate finite difference approximations 
for the entire I-beam. 

We do this by plotting the amplitudes versus the corresponding values of 
k¢for the fundamental, first, and second modes. The presence of any peaks 
would indicate a resonant type behavior. If any of these peaks corresponds to 
a particular value of ky for a clamped or fixed plate we say that for that value 


of kg (and hence @), the I-beam behaves in a manner similar to a plate with 
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clamped or fixed boundary conditions. This gives us a possible range for ky 
values on which to concentrate when looking for resonant frequencies. 

The values of ky are given in Table 2. The quantities in parentheses are 
the modes that propagated for a particular value of ky which is determined by 


the maximum integer value n can take such that By, is real (Ref. Equation 71). 


TABLE 2. kg VALUES 


— SYSTEM 
CATEGORY | ORIGINAL 


1 (0) 1.2516 (0) 0.5265 (0) 
4 (0, +1) 3.3446 (0, +1) 2.026 (0) 
7 (0, +1, +2) 6.5572 (0,+1,+2) | 4.558 (0, +1) 





We divide our values of ky into three categories (for identification 
purposes). In the first category we include the values of ky corresponding to 
the first resonant frequency of the clamped and fixed plates and the values of 
kg for our original system which allows only the fundamental mode to 
propagate. In the second category are values of ky that correspond to the 
second resonant frequency of the clamped and fixed plates as well as the value 
of kr which allow only the fundamental and first modes to propagate. The 
third category corresponds to the third resonance of the clamped and fixed 
plates and the first three modes in our original system. Also included in this 
final category are the values of ky = 7.5, 8 and 9 which will allow us to study 


the behavior of the second mode at higher frequencies in greater detail. 
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In Figure 19 we show a plot of ky values versus the amplitude of the 
fundamental mode. Computed values are shown in black. A cubic spline of 
the points involved is shown in red. We must note however that this spline 
is only a possible representation of the behavior of the amplitude for different 
values of kf, since time constraints did not allow us to conduct a more 
comprehensive set of simulations from which we could obtain an accurate 
picture of the behavior. 

It can be seen that for k¢ = 0.5065 and 2.026 (the points labeled FF1 and FF2, 
the I-beam is exhibiting a resonant type behavior. These values correspond to 
the first and second resonant frequencies of a fixed plate. It can also be seen 
that for values of ky greater than 3.5 there is little variation in the amplitude 
of the fundamental model since we are now past the first cutoff value (after 
which three modes propagate) of ky= 1. This leads us to believe that at these 
higher frequencies most of the activity takes place in the higher modes. 
However the total energy calculation shifts from fundamental to first and 
back to fundamental as can be seen in Table 3. 

In Figure 20 the only resonant behavior is exhibited for the value 
kt = 3.446 (the point labelled CC2). This frequency is just past the first cutoff 
value (z) and the amplitude of the first mode is double that of the 
fundamental (compared with the corresponding value in Figure 19) which 
confirms that at higher frequencies energy is being propagated in the higher 
modes. 

In Figure 21 we plot the amplitude of the second mode against ky. The 
only resonant type behavior which exists here is for the ky = 8, but since the 


amplitudes involved are so small it would be difficult to draw an accurate 
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conclusion about the existence of a resonant frequency. Another aspect we 
investigate is the proportion of the total energy carried by each propagating 


mode which is given by the formula 


Z 
pM 


M 
4h 
p=-M 


(159) 


where E,, is the energy of the nth mode and A, the amplitude of the nth mode 


(Kinsler, 1982, p. 110). These quantities are summarized in Table 3. 


TABLE 3. ENERGY CONSIDERATIONS 


G7 2 82.43 
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We can see that the proportion of energy carried by the fundamental mode 
varies from 100% to 12.04% as kg goes from 2.026 to 3.3446. This can be 
accounted for as a redistribution of energy that takes place as the frequency of 
the waves passing through the fluid wave guide pass the first cutoff value of 
k¢= 2. Similar arguments hold for higher values of ky. 

Finally we take a brief look at the shear strain field generated in the solid. 
We do this for two reasons, one to check the behavior at the corner nodes and 
two, to find out where the maximum shear strain occurs. 

When treating the corner points we neglected applying the stress free 
boundary conditions there assuming this would have no adverse effect on 
the behavior of the solid. Were this assumption invalid we would expect to 
observe singularities or other irregular behavior. Provided in Figures 22, 23 
and 24 are snapshots of the solid at different time levels which display no 
unusual behavior at the corner nodes, leading us to accept the assumption as 
valid. 

From an engineering standpoint we are interested in where the 
maximum shear strain occurs for possible failure analysis. As can be seen 
from Figure 25 the maximum occurs along the transverse borders of the 
I-beam cavity. We must note that the amplitudes and frequencies of the 


acoustic pressure waves involved in our study are not comparable to those 
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which would cause permanent deformation or failure, but which could result 


in fatigue cracking if subjected to prolonged periods of stress. 
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V. CONCLUSIONS 


Several problems arose in implementing the numerical scheme. It was 
found that the non-local radiation boundary condition did not completely 
annihilate the initial wavefront and a small amount of reflection took place, 
but over time the effect of this reflection became negligible. Evanescent 
modes did not decay sufficiently and were reflected at the boundary adding to 
the overall noise of the problem. 

When calculating the amplitudes of the propagating modes at steady state 


we applied the formula 


1 
A= | p yaneay (161) 
i 


which should yield consistent results for any y values in the fluid domain. In 
the neighborhood of the fluid/solid interface this was not the case as can be 
seen from Figure 26. We believe this is due to the presence of high frequency 
components in this region and having too coarse a mesh to effectively 
evaluate the integral there. 

With a stepsize h of 1/40 our truncation error was on the order of 107°. 
To increase the accuracy we can perform Romberg extrapolation or decrease 
the stepsize, thereby increasing the dimensions of the matrices involved. The 
latter was not an option due to the construction of the code and machine 
limitations, in that a simulation which involved 10 to 15 thousand time steps 


took from 10 to 12 hours to perform. Increasing the size of the matrices 


iS. 


involved, thereby increasing the required time, would not make it a suitable 


code for experimentation and timely results. 


Amplitude of fundamental mode 
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Figure 26. Amplitude of Propagating Modes in the Y-domain 


Aspects of the problem which deserve further study would be the use of 
obliquely incident waves vice normal as was used here. Variation of the 
cavity size and its effect on the amplitude of the propagating modes should 
also be considered. As a simplifying assumption the cavities were considered 
to be void. It would be realistic to expect that they may contain fluid 


(seawater) or an acoustic dampening material. Modifying the model to 
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account for such conditions warrants further study. A final aspect that could 
be considered is a resonance analysis. 

Throughout all of our finite difference approximations we were able to 
maintain second order accuracy in space and time. We did not have to resort 
to the use of pseudonodes outside of our fluid-solid domain. Finally the 
possibility of being able to establish an acoustic signature for a double-hulled 
structure could open up a new avenue of submarine detection for which this 


thesis could be considered a starting point. 


is, 


APPENDIX A. VON NEUMANN STABILITY ANALYSIS FOR THE 2-D 


SCALAR WAVE EQUATION 


The general form of the scalar wave equation being used is 


ate 


O7u atu ou 
2 ax ay? 


whose equivalent difference equation, given Ax = Ay =h is 


( n+1 n rt) 1 
a i 
cf At? 1,] 1,} 1,] 


The error function takes the form 


u = EE ep h 


e-T) 


1 
n n A —lur . 
7 =o (Mle, uf +uf, j)+ h? (uj 2a) +2} 


(A-2) 


(A-3) 


This is substituted into Equation A-2 and common terms are cancelled to give 


oe —2+ oad = alle ~2+67Hh) 5 (ef -2+e 7%) 


Using the following identity, 


10x +e tax 
COS OX = 


letting Bh = yh and defining p = cAt/h, Equation A-4 reduces to 


i) 2 p*(4cos Bh — 4). 


Multiplying Equation A-6 through by € and collecting like terms gives 
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(A-4) 


(A-5) 


(A-6) 


Ee ~ 2€(1-2p?(cosBp-1)}+1 =0 


or equivalently 
ae 2e(1 Se sin?) 1=0. 


The roots of this quadratic are 


= 
E182 = 1- 4p? sin®( 2° i{1- 4p? sin?) Fis 


For stability we require 


E<] 


which forces 
ph\\ 
{1 ae sin’) ~4<0. 


Solving for p* in Equation A-10 gives 


» 1 
, o——— 
asin?( E) 
w 


which reduces to our final stability criterion of 


Or 


(A-7) 


(A-8) 


(A-9) 


(A-10) 


(A-11) 


(A-12) 


APPENDIX B. VON NEUMANN STABILITY ANALYSIS FOR THE ELASTIC 
WAVE EQUATION 


In this section we give a brief outline of the stability analysis required for 


the elastic wave equation, whose general form is 





20 0 
gu, 204 Bs ae B-1 
he + cf ay? +(cf cf | andy OF (B-1) 
dv 90° 9\ a-u — a7v 
oo =) = = Bae 
Ta cL ay? ( L cf | axdy WF (B-2) 


We use the following error functions 
u= Ur ellBp+ mi) and v= ver ellbptm)y (B-4) 


which are substituted into Equation B-3, common terms are cancelled and 


complex exponentials are gathered into trigonometric quantities to give 


—4r *[<fsin? Bn sin? — it u- r (cf - cf \(sin Bhsin yh)V = (€-2+87)-U(B-5) 


Following the same procedure for Equation B-2 we obtain 
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—4r (cf sin? e +C7 e sin? P\y- r (cf - ch \(sin Bhsin ym)u =i —2+ E-!).V.(B-6) 
Equations B-5 and B-6 can be written in matrix form which is 


aii (2 ae +c&sin? # -r? (cf = cf )sin Bh sin 1 at a 
kv A 


ae (cf - sina sin ¥1 —4r iG sin? — + Cf 2 sin? tie 


where A = &€-2+ 1. The eigenvalues of matrix in Equation B-7 are 


Os eS gate 1 
ICT Cr = 
Pie (cf +¢F (cos Bh + cos Yi —2)+ sana (a — yh) +cos(Bh + v1) - 2\° ‘ 
(B-8). 
they take on a maximum value when fh = yi = a, and we obtain 
Ay =Ao = —4r? (cf + cf (B-9) 
We are now left with the identity 
See ape. =-4r7 (cf + cF (B-10) 
or 
E2 + (4°? (ci +cf)- 2)E+1=0 . (B-11) 


For stability we require 


which forces 
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OF 


(0 


2 


2 
ct + cf) - 2) -~4<0 





APPENDIX C. FINITE DIFFERENCE FORMULAE FOR THE EQUATION OF 


MOTION AND BOUNDARY CONDITIONS 


This appendix contains the finite difference approximations of all the 


equations required for the simulation. 


Solid 
1) Elastic Wave Equation 


2 2 
n+1 At n n At ( n n | 
Us. =—mo—aelUs, tus, -|t+—aamlue sy tu 
if kph? Neel ie} ken? loiard id 


1 1 n n Ks : 
“ 4h? ES , Jj ~ Oi j+1 Ci, j-1* ja) 


= (v! +01 Ae vr. 4 tue 
ey keh? feeily| ley keh? Ljjerll el 





2 
At 1 1 ( n n n n | 
+— ~ Uta tragic eT a 
{a 7 Celie ijl Heli tal Te fae 


81 


(Gey 


(C - 2) 


Stress-free Boundary Conditions 


Surfaces with f = G | 

ie Gara ar)) Ci eon ena) 
Site = + (et, j=l Vis, ia} = la Z “ae 7 oF) al. 
(Se eg Bree et agin 


At?7{1 1 At#?{1 3 
Satara Cone —uP je] - ae ae q uF j-1) (Can 


—] 
Surfaces with nf =(( 
2 2 2 
n+1 = _ 2At 1 1 n an n—] 2At n At ( n n 
Ui -[2 om fare Ui Uj eR Mist,j * yy Us ina t Up iy 
ps 2 
At? Ta ( Pa ‘ Ai lo ( “ * 
—-—~| ~-— |v... ,-0:,,.; +—~| ~-— |v... -0; (C—5) 
a(S 7 tee fod ictal etal oh2 ee k2 Ti Lfao 
2 2 2 
fer _ 2At as ey n mer) 2At” in At (o/ n 
vi -[2 ao Fare i ij vi Pe Uist ap U5 a1 + Uj jy 


2 2 
AIS) Ae ee Ary st (us Lye = 
a2 F ke lta ja ee) he Ee | a tla (C—6) 
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0 
Surfaces with nf = | 


2 2 
2At?( 1 1 fin ce AVE At 
yntl Boas tao | (4. i) ey,” +a (u -+ur ’ 
i} h2 k? ké a 1,] h?ke jal hk aly ih] 
46a 1 (on | ot Ea 13 \on on, ) 25, 
oh2 k? ka Paleo i+1,)-1 rh k? k2 i+1,j lb 
ynt! | 2- 2At" a Us _ pl, 240 Us 4 (v +u! 
1,J Pe kf ke 1,] 1,J hk Pia hk itl] Jelly 


2 Z 
wad) 1) , P Aid| clon : 
+S eae (eae) Se ea [eH a 


0 
Surfaces with ni = a 


2 2 2 
oe | 2 — ——____| —_ + —__ | ly". —y. +a; i: +. 
1,] | h2 = k2 i ME) hk? amas nk? Pee tL 


- 2 
ee es Cee ; 
ae Yi raya) 3 eho, Ui-1,j (C-9) 


antl _ es 11 yn ofl Zoe ‘ AP fon wn 
1,J 7 h? ke ae i hk Cel h?ké Ui obi Ui4] a 


2 2 
ia 1 7A) Gn. | (aes aes. 
-25(5- act = uh jaa) ety : 5 uta, -4h i” 
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Corners 2 and 4 


2 
n+1 2At°} 1 1 n n-1, At ( n n |+ At (ut n | 
Uo. S12 — +— {ju.--u:. +—— I. ty: —U.. 
Ns) 1 ars k2 1,] 1,] h2k? fools, lh) keh? i,j+1 ~“1,j-1 
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Corner nodes 1, and 3 
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At 1 1 nN fl n n 
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At 1 1 n n n n 
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Fluid 
1) Wave Equation 


ntl _|>_ 440° ra pte, AL (yn tot. ay! 4 lt he=15) 
Pi = Kh Pig Pi,j Kh Menibe ele reer i, jk 


Radiation Boundary Condition 


—] 
n+1 _ At 2At* n ie a At n-1 
a [8 4 Es Pi, ja TP, t sep a 


AAt? At? At? , 
21.2 21.2 21.2 
kth kh kfh 
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0 
At? 
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kph* 
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k¢h k¢h kth 


85 


0P JP oP oP 


APPENDIX D. COMPUTER CODE 


function basel - base for the ‘I - beam ’ shaped domain 
written by: Lt. Hugh Mc Bride USNR 

date : Apr 92 

constructs a grid in the shape of an I beam 


function [b,rows, pcoll,pcolr,fill,coll,colr] = basel(n) 


FP AP dP dO HO OP OP OP OF OP OP Of AP AO OP OP 


de oO 


variables 

rows - identifies the rows to be zeroed to form the I beam 

coll - identifies the columns to be zeroed to the left of the center spar 
colr - identifies the columns to be zeroed to the right of the center spar 
pcolr - identifies the columns to be zeroed to the right of the center spar 
including the boundary values. 

pcoll - identifies the columns to be zeroed to the left of the center spar 
including the boundary values. 

bb - building block of the correct size 

cc - dimension of the blocks to be zeroed out on either side of the center 
spar 

S - variable to adjust the size of the zero blocks and keep it symmetric 
fill - block of zeros of appropiate size to fill in the spaces on 

either side of center spar i.e. zeroing out,at every time step. 

m- no. of divisions in the half length of the domain. 

cr - variable used to pick the columns to the right of the spar 


build a grid of the correct size. 
bb = ones(2*n+1); 


the basic variables required to build the I - beam shaped grid 
m= ntl; 
S = log2(n); 
cc = (s-1)*(log2(n)); 


determine the row numbers to be zeroed out. 
rows= m-(cc-1) :m+(cc-1) ; 


columns to the left of the center spar, pcoll includes the boundary points 
coll does not. 

coll = 1l:cctl1; 

pcoll = l:cct2; 


columns to the right of the center spar, pcolr includes the boundary points 
colr does not. 


cr = 2*n+1 - (cc); 
colr = cr:2*n+1; 
peolr = cr-1:2*n+1; 


construct the block of zeros used to fill in the spaces on either side 
of the center spar at every time step. 

e = size(rows); e = e(2); 

@ = size(coll); d = d(2); 

fill = zeros(e,cctl); 


86 


bb(rows,coll) 
bb(rows,colr) 


b = bb; 


5 oan 
fill; 
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function bndfnx - boundary facing in the negative x direction 
i.e unit normal = (-1,0) 

written by: Lt. Hugh Mc Bride USNR 

date > Apr 92 


bndfnx applies the traction free boundary and the elastic 

wave equation along the boundary of the I-beam with unit normal 
(-1,0) the technique used is that as developed by Ilan and Lowenthal 
and is not discussed here. 

the columns containg the nodes on the surface in question and 
the required neighbours (those one column in from the surface) 
are picked off from the current and old values of u and v 

the shifted and weighted with constants from the vectors 

cunx (Constants for U-values for the Negative X boundary) 

and cvnx (Constants for V-values for the Negative X boundary) 
and inserted in the correct position of the updated u and v. 


PO AO AO AO HO AO AO AO AO HAO HO HO HO OP HO OO 


function [un ,vn) = bndfnx(uc,vc,u0,vo,un,Vvn,rows,pcoll,c,d); 


Variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values of v 

rows : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 
the nodes on the boundary. 

pcoll carries out the same function as rows for 

the columns, and the contains the column location of the 

of the boundary facing the negative x direction 


PP AO HO AO ' AO HO HO AO OO 


cunx; 
Cvnx; 


a oo 
m%O 
t ot 


(13,33) = size(rows); 
mrows = [rows(13)-1 rows rows(j3)+1]; 
[14,34] = size(mrows); 
[15,35] = size(peoll) > 


cul cu2 cvl cv2 co cov contain the necessary u and v values 
for our calculations 


of oo 


cu2 = uc(mrows,pcoll(j5)+1); 
cul = uc(mrows,pcoll(j5)); 


Cv2 = vc(mrows,pcoll(j}5)+1); 
CVI = Ve(mpows, pcoll (75); 


co = uo(mrows,pcoll(j5)) 


; 
Cov = vo(mrows,pcoll(j5)); 


% the updated values are calculated 


uci = c(1l)*cul(2374-1) - co(2:34=1) + ¢(2)*cu2z2- 42)... 


+e(3) *(eul(3?34)04 Cull? )4—=2)5 
+0 (4) *(cv2(1: 34-2) g= Cvg@i- ag)... 
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+O(5) (evils: ya) — CV1(1: 34-2) ); 


Mete— d( 1) *evl(2:34=-1) — cov(2:34=-1) + a@(2) *cv2(2: 34-1)... 
Tas) (ev l( sega) + CV1(1°34=2)) 22. 
+d(4)*(cu2(1:34-2) - cu2(3:34))... 
+€(5)*(cul(3:34) - cul(1:34-2)); 


% and put in their proper place in un and vn 
pb= size(mrows) ; 


pbc = mrows(2:pb(2)-1); 


Wmeebe,pcoll(j35)) = ucl; 
vn(pbc,pcoll(j35)) = vel; 
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function bndfny - boundary facing in the negative y direction 
i.e unit normal = (0,-1) 

written by: Lt. Hugh Mc Bride USNR 

date : Apr 92 


bndfny applies the traction free boundary and the elastic 

wave equation along the boundary of the I-beam with unit normal 
(0,-1) the technique used is that as developed by Ilan and Lowenthal 
and is not discussed here. 

the rows containg the nodes on the surface in question and 

the required neighbours (those one row in from the surface) 

are picked off from the current and old values of u and v 

the shifted and weighted with constants from the vectors 

cuny (Constants for U-values for the Negative Y boundary) 

and cvny (Constants for V-values for the Negative Y boundary) 
and inserted in the correct position of the updated u and v. 


GO GO dO dP OP dO dP OP OP OO BO dO OP CO OO OO 


function [un ,vn] = bndfny(uc,vc,uo,Vvo,un, vn, rows, pCO! pcoln, ca 


variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values of v 

rows : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 

the nodes on the boundary. 

pcoll and pclor are both required as there are two regions,one on 
either side of the center spar of the I -beam which require 

our attention and they contain the location of the nodes in question 


OO AP WP HP AO CO OO AO GO CO 


cuny; 
cvny; 


oe oo 
on 
tf il 


(11,51) 
(12,j2] 
(Sr se] 


$1ze(pcoll); 
size(rows) ; 
size(uc); 


rl** ru* rv* and ro* pick off the rows on either side of the 
center spar of the necessary u and v values. 


of oO 


rlul = uc(rows(j2)+1,pcoll); 
= uc(rows(j2)+2,pcoll) ; 
ru2 = uc(rows(j2)+2,pcolr) ; 


Eu uc (rows (j2)+1,pcolr); 
rlv1l = vc(rows(j2)+1,pcoll); 
rlv2 = ve(rows(j2)+2,pcoll); 
rv2 = ve(rows(j2)+2,pcolr) ; 
rvl = ve(rows(j2)+1,pcolr); 


rol = uo(rows(j2)+1,pcoll); 
ro = uo(rows(j2)+1,pcolr); 

rolv = Vol rows(72)4 1 peel, 
rov = VelTovs( 72) +1, peclr)- 
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% the updated values are calculated 


ie e(l)*rlul(2:jl—1) = rol(2:jl-1) + ¢(2)*rlu2(2:j1-1)... 
Peseta (321) + rlui(1:31-2))... 
Tea pe(rive( isa l=2) — riv2(3:31))... 
+¢(5)*(rlv1(3:jJ1) - rlvi(i:j1-2)); 


esi (23 9l—2j— ro(2:j1l-1) + c(2)*ru2(2:j1-1)... 
money (ral(3sj1) + rul(ivjl—2y)..: 
mea) * (V2 (1. a2) - rv2(3:31))... 
Hema) <(rvid(ougioe — FVvVl(1:j)1=2)); 


vl = d(1)*rlv1(2:31-1) - rolv(2:31-1) + d(2)*rlv2(2:j1-1)... 
meek (Civiltssj1) +. .rlv1(1:j1-2))... 
mene (rlu2 (il. ji=2) = rlu2(3:31))... 
fa5)*(riul(3:j1) - rlul(1:jl1-2)); 


ieeai( i) *rvl(2:j1-1) - rov(2:ji-1) + d(2)*rv2(2:j1-1)... 
memes ye (rV1(3:31) + rvl(i:j)1-2))... 
ora) *(ru2(1l:jl-2) = mu2(3:j1))... 
ta(5)*(rul(3:j31) = rul(i:j31-2)); 


% and put in their proper place in un and vn 
pb= size(pcoll)j; 


pbl peoll(2:pb (2) -1); 
pbr meoln(2:pb(2)-1); 


to tt 


un(rows(j2)+1,pbl) = ul; 
un(rows(j2)+1,pbr) = ur; 


vn(rows(j32)+1,pbl) = vl; 
vn(rows(j2)+1,pbr) = vr; 


% the same procedure is reppeated for the ‘bottom’ of the I-beam 


Mammo Sc—-1) = c(l)*uc(1,2:sc-1) = uo({1l,2:sc-l1)... 
mene) *UC(2,2:Sc-1)... 
pees *(UC(1,3:sc) + uc(1,l:sce-2))... 
Pea yC( 2, 1:Sc=2) = ve(2,3:sc))... 
fec>)*(vc(1,3:sc) = ve(1,1:sc-2)); 


emeeeesc-1) — G(1)*vc(1,2:sce-1) - vo(1,2:sc-1)... 
tmar2)*vc(2,2:sc-l)... 
momo (VC(1,3:sc) + ve(1,1:sc-2))... 
momar *(uc(2,1:se-2) =- uc(2,3:sc))... 
mos) = (C(1, 3:sc) = uc(1,1:sc-2)); 
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% function bndfpx - boundary facing in the positive x direction 
$ i.e unit normal = (1,0) 
%$ written by: Lt. Hugh Mc Bride USNR 
+ date [Api oe 
% 
% bndfnx applies the traction free boundary and the elastic 
% wave equation along the boundary of the I~-beam with unit normal 
¢ (0,1) the technique used is that as developed by Ilan and Lowenthal 
% and is not discussed here. 
% the columns containg the nodes on the surface in question and 
% the required neighbours (those one column in from the surface) 
% are picked off from the current and old values of u and v 
%* the shifted and weighted with constants from the vectors 
% cupx (Constants for U-values for the Positive X boundary) 
%* and cvpx (Constants for V-values for the Positive X boundary) 
%* and inserted in the correct position of the updated u and v. 
function [un ,vn) = bndfpx(uc,vc,u0,Vve,un, vn, cows, peolr,c ca) 
% variables 
% un : updated values of u vn : updated values of v 
%* uc : current values of u ve : current values of v 
% uo : Old values of u vo : old values of v 
% rows : used to identify the elements of the matrix which 
% are zero for all times they also contain the row location of 
* the nodes on the boundary. 
% pcolr carries out the same function as rows for 
%¢ the columns, and the contains the column location of the 
% of the boundary facing the positive x direction 
% c = Cupx; 
% d = cvpx; 
{i3,}3) = size(rows); 
mrows = {rows(i3)-1 rows rows(j3)+1); 
{i4, 34] = size(mrows); 


$ cul cu2 cvl cv2 co cov contain the necessary u and v values 
* for our calculations 


cud = uc(mibows, pecln(2)—1); 
cul = uc(mrows, peolr(1))- 
cv2 = vc(mrows,pcolr(1)-1); 
Cvl = vce(mrows,pcolr(l1)); 


co = uo(mrows,pcolr(l1)); 
cov = vo(mrows,pcolr(1)); 


° 


%* the updated values are calculated 


ucr = 6(2) *eu1(2:)45=1)h —.60( 2): 4-1)" + 6( 2) *eu2 (2: 94-1)) 
+@( 392 (Cul (3 34) ten) (4) 44 —2 ee 
+6(4) * (CV2 (eg 4—2 )a— eve (354) ) ee 
+0(5)*(CGV1(3: 94) -="evl( ts 4—2))); 


a2 


weeme—  1) *CV1(2:954-1) = cov(2:34-1) + d@@)*ew2 (2634-1) cm. 
Fos) = (evils: 44) + cvl(1:34-2))... 
Fa (4)*(cu2(1:74-2) = cu2(3:34))... 
+4(5)*(cul(3:34) - cul(1:34-2)); 


% and put in their proper place in un and vn 


pb= size(mrows) ; 
pbc = mrows(2:pb(2)-1); 


un(pbc,pcolr(1)) = ucr; 
Virgepc,pcolr(1)) = vcr; 
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function bndfny - boundary facing in the positive y direction 
1.e unit normal = (OF) 

written by: Lt. Hugh Mc Bride USNR 

date : Apr S2 


bndfpy applies the traction free boundary and the elastic 

wave equation along the boundary of the I-beam with unit normal 
(0,1) the technique used is that as developed by Ilan and Lowenthal 
and is not discussed here. 

the rows containg the nodes on the surface in question and 

the required neighbours (those one row in from the surface) 

are picked off from the current and old values of u and v 

the shifted and weighted with constants from the vectors 

cupy (Constants for U-values for the Positive Y boundary) 

and cvpy (Constants for V-values for the Positive Y boundary) 
and inserted in the correct position of the updated u and v. 


function (un ,vn) = bndfpy(uc,vc,u0,vo,uUn,Vn,LoWws, pCO!!! ,pco!'., cyan 


BP NP AP AO OO DP AO OP CW OP 


% 
% 


Cath a) 
(12,2) 


rel} biel 
EluZ 
ruz 
rul 


yal 
riv2 
EVZ 
5 a 28 


variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values ct vy 

rows : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 

the nodes on the boundary. 

pcoll and pclor are both required as there are two regions,one on 
either side of the center spar of the I -beam which require 

our attention and they contain the location of the nodes in question 


size(pcoll); 
slze(rows) ; 


ft Ul 


r]** ru* rv* and ro* pick off the rows on either side of the 
center spar of the necessary u and v values. 


uc (rows(1)=—1,pcoll); 
uc(rows(1)<-2,pcoll); 
uc( vews (1) =27 peolm); 
uUC(TOWS (1) -1, peol,):- 


ft Wl 


ve (rows (1) <1, peel!) > 
ve(rows(1)<2,pcoll) ; 
ve(rows(1)-2,pcolr); 
Ve(rows (1)=-1,peolr); 


rol = uo(rows(1)-1,peoll); 
ro = uo(rows(1)-1,pcelr); 

rolv = vo(rows(1)<-1,pcoll); 
rov = vo(rows(1)-1,pcolr) ; 
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%¢ the updated values are calculated 


femeeec (1) *r1ul(2:j31-1) - rol(2:31-1) + c(2)*rlu2(2:31-1)... 


Peto yx(rlul(3:9l) + riai(1:31-2))... 
Terayed(rilv2(l2jil—2) = rlv2(3:j31))... 
ese (rivi(ss gl) = rivi (1:j31-2)); 


meee. 1) *rui(2:j1-1) = ro(2:j31=-1) + c(2)*ru2(2:j31-1)... 
memes (rul(3: 71) + rul(1:j31-2))... 
memayet(rv2(1% )1=2) = rv2(3:3)1))... 
+¢0(5)*(rv1(3:j1) - rv1(1:31-2)); 


Seemed bh) xr lvl(2:j1-1) = rolv(2:31-1) + d(2)*rlv2(2:j31-1)... 


Tomsk (rivl(3:jJ1) + rlvil(1l:31-2))... 
+d (4) *(rlu2(1:31-2) - rlu2(3:31))... 
met) *(rlul(3:31)- rludl(i:31=2)); 


Seal ( 1) *®rvi(2:j1<-1) - rov(2:31-1) + d@(2) *rv2(2:j1-1)... 
+d(3)*(rv1(3:j1) + rv1(1:31-2))... 
moma) *(ru2(1i:j1-2) = ru2(3:j1)).. 
Bares) *(rul(3:j1) = rui(is¢j1-2)); 


oe 


and put in their proper place in un and vn 
pb= size(pcoll); 


Ppt = pcoll(2:pb(2)-1); 
=mmeColr(2:pbp(2)-1) ; 


Mme rOwS(1)-1,pbl) = ul; 
timenows(1)-1,pbr) = ur; 


wn (rows(1)-1,pbl) = vl; 
vn(rows(1)-1,pbr) = vr; 


D5 


function cmodck - checks the amplitude of the propagating mode 
written by: Lt. Hugh Mc Bride USNR 
date ; Apr 92 


cmodck is the driver program for the problem.All the vaue of the constants 
are defined here and all quantities are scaled before they are fed into the 


of oP GP WO AP CO 


decks are cleared before calculation 
record(’erase’) 
clear 
clg 
axis(‘auto’) 


oe 


ax = input(’step size = ‘); 
a= input(’ length scaling factor = ‘); 
input(’ time scaling factor (omega) = ‘’); 


O 
J 
iQ 
il 


constants 

ct: tranverse velocity of the solid (steel) 
cl: longitudinal velocity of the solid (steel) 
dnss: density of the solid (steel) 

cf : speed of sound in fluid (seawater) 

dnsf : density of fluid (seawater) 

epss: ratio ot fluid to solid densities 


GP AO WO DO oP oP OO 


ct = 3200; cl = 5900; cf = 1500; dansf = 1026; dnss = 7700; 
epss = dnsf/dnss; 


% the determination of the scaled variables 
% axs scaled distance 

% ats scaled time 

axs = ax/a; 

kf = omg/cf;kt = (omg*a)/ct; kl = (omg*a)/cl; 


% datsl and dts2 are the stability criteria for the solid and fluid 
% always choose the minimum 

@tsl = (Ki*axs), (squu(l + (ct-2/cl-2)) 5 

ats2 = .5%*((kKf*a*daxs) /sqrt(2)); 


[dtsl dts2] 


ats = input(’ desired time step = ’); 
nn = input(’ number of timesteps = ’ ); 


% st - when to stop building the radiation boundary condition 


for in the construction of the matrix A each loop thru the construction 
% cylce allows another mode to propagate 


st = stop(kf); 
1 = sqrt(-1); 
 X a vector the length of the domain used in several places 
% i.e. when integrating etc. 
%$ mm ~ the mode being checked O-fundamental etc. 
* nt - no of intervals in domain 
% tt weighting factor for the trapezoidal rule 
X = =Tsaxs:]. ; 
m= size(x) ; m= m(2) ; ml = zeros(m) ; 
mm = 
nt = (1/daxs); 
EC) ="(m—=P ye © 
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{un,vn,mn,zZ0O} = cwv4(m1,kl,kt,kf,dts,dxs,x,nn,st,epss,tt,nt,mm); 


% 
% 
% 
% 


un 
vn 
mn 
zO0 


displacment in the x dirn 
displacment in the y dirn 
displacment of the fluid 
vector containing the amplitude of the propagating mode 


7 


% function coupl - couples the two media at the fluid solid interface 
% written by: Lt. Hugh Mc Bride USNR 
% date : Apr 92 
% 
% coupl applies the traction boundary conditions at the fluid solid 
% interface where the shear component is zero but the normal component is 
% given by tau(yy) = -p(total) ,which causes the I beam to deform 
% and causing the propagation of scattered pressure waves to propagate in 
% fluid. To prevent cavitation at the interface we apply the inviscid 
% form of the Navier-Stokes equation ( or Euler’s equation ) at the 
% the boundary 
function [un ;vn,mnpd]}] = coupl(uc,vc,uo,vo,un,Vvn,c,d,k,dts,dxs,all, biligepecum 
% variables 
% un : updated values of u vn : updated values of v 
t uc : current values of U ve : current values of v 
* uo : old values of u vo : old values of v 
%$ ¢ = cupy; 
* d = cvpy; 
1 = sgrt(-1); 
[sr sc] = size(uc); 


% the forcing function 
time = exp(-i*k*dts) *ones(1,Ssc) ; 


¢ the u component requires no modification and is treated in the 
% usual fashion. 
un(sr,2:sc-1). = c(1)*uc(sH,2:Sc-1) =sne(sr, 2: Sc= ieee. 
+ C(2)iFUC(SE-l, 22Se— 1)... 
+e (3) (We(sr, 3 7Se)r +uc(Ssr, 1 Se-2)) 4. 
+ C( 4) * (Vc (Sr-)7 1 3Se22)— Vo(Sre— 13 Se) jn. 
+ C(5) *(VeCl(Sr,3 See -— VeE(S., 1 .Sse-2). 


% the periodic boundary condition for u 
un(Sr,1).2=-c(1)4ue(sr, 1) = Ju0(sr,1)... 
+ CC2)*uc (Sri) lua. 
+ ¢C(3)*(e(Sr, 2) 3+ uc (Sr,Se—-l ra 


+ c(4) *(ve(Sr-1,Ssc-1) - ve(Sr-1,2))... 
+ ¢C(5)AGve(sSr,2) >> VeE(sr,se-1))). 
un(sr,Se)"= un(sr;l))-; 


the normal component of stress plus the incident,reflected 
% and scattered pressures 
vn(sr,2:scr-1) = d(1)*vc(sr,2:sc-1) - vo(Sr,2:sc-1)... 
+ d(2)*vc (sr=l2+sc-1).... 
+d:(3) *(Ve(Sr, 3 7ser + VE(Sr, 1 -Se-2)))... 
+4) = (UC (Sr—1,.1:Se—-2)) — ue(Si—-ie 3:Se))....- 
+0(5) *(uc(sr,3:sc) = uc(Srise-2))... 
+ ((2*dts*2) /dxs) *(2*epss*time(2:sc-1) + epss*al1l(2:sc-1)); 


% the periodic boundary condition for v 
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fest ,.) = d(1)*ve(sr,1) - vo(sr,1)... 
fmcleeye VC (Sr=1,1)... 
Home *(VC(Ssr,2) + ve(sr,Ssc-1))... 
emtta) = (UC(Sr—1,Se=-1) = uc(Sr-1,2))... 
Faro) *(uc(sr,2) — uc(sr,sc=-1))... 
+ ((2*dts*2) /dxs)*(2*epss*time(1) + epss*al1(1,1)); 


wn(sr,sc) = vn(sr,1); 


% the compatability condition. 
moeee=— (—-(2*axs) /(dts“2))*(vn(sr,:)- 2%ve(sr,:) +vo(sr,:)) + b11(1,:); 


oe 
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function cwv4 - coupled waves version 4 

written by: Lt. Hugh Mc Bride USNR 

date > Apr 92 

cwv4 is the main program which couples the behaviour of the ‘/I-beam’ shaped 
solid medium with the fluid medium. The elastic wave equation and the boundary 
conditions ,periodic and traction free are are satisfied for the solid, while 


the scalar wave equation and the periodic and radiating boundary conditions ar 
applied to the fluid. 


The general steps of the program are as follows 


1. The basic parameters are determined ,that is the size of the domains 
as passed to it by the driver program cmodck.m 


2. All the global variables are calculated for both the fluid and the solid 
including the matrices required for the radiation boundary condition. 


3. vl and v3 are vectors used to find the amplitude, vl is the 
weighting vector 1/2 111.... 11/2 for the trapezoidal rule 
and v3 = exp(i*n*pi*x) , the othhogonal vector. 


4. The I beam is built by basel 


5. The initial values of the u and v for the solid (un,vn,uc,vc,uo,vo) 
and m (mn me and mold) for the fluid are set to zero. 


6. The elastic wave equation and the periodic boundary conditions 
are applied to u and v 


7. The boundary conditions for the solid are applied 
8. The fluid and solid are coupled 
9. The freespace portions of the I beam are zeroed out 


10. The scalar wave equation and the periodic boundary conditions for 
the fluid are applied. 


11. The values of the amplitude are calculated and accumulated 


12. u, v and m are updated ,that is the new values become 
the current values and the current values become the old values. 


function [(v,un,vn,mn,Z2] = cwv4(m,kl,kt,kf,dts,dxs,x,nn,st,epss,tt, nem 
variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : Old values of u vo : old values of v 

rowS : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 

the nodes on the boundary. 

pcoll and pclor are both required as there are two regions,one on 
either side of the center spar of the I -beam which require 

our attention and they contain the location of the nodes in question 
kl - scaled longitudinal speed 
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kt - scaled transverse speed 

axs - scaled spacing 

dts - scaled time step 

x - vector of length of the interval 

nn - no of timesteps 

st - stopping criteria for radiation boundary condition 
epss - ratio of solid density to fluid density 

nt - no of intervals in domain 

mm - prpoagating mode under investigation ,mm = 0 fundamental 
mm = 1 first and so forth 

z? - amplitude of propagating mode 


JO AP AP AO AP OP AO AP AP AP HO 


Step 1 : r,c the dimensions of the domain, r-1 ,c-1 and n-l 
the dimensions of the interior 


Je do 


[r c}] = size(m); 
Nem-erel; ©£r = r-l;¢c = ¢c-1; 
axis (’xy*) 


% Step 2: The global variables 


% For the solid 


mao = rhov(kl,kt,dadxs,dats); ces8= (2° — 2*rho(1)-2*rho(2)); 
eunx = ucnx(kl,kt,dadxs,dts); cvnx = venx(kl,kt,daxs,dats); 
eumy)— Ucny(kl,kt,dxs,dts); cyny = veny(kl1,kt,daxs,dts) ; 
eux — UCDX(K1,kt,dxs,dts); cvpx = vepx(kl,kt,dxs,dts) ; 
cupy = ucpy(kl,kt,dxs,dts); cvpy = vcepy(kl,kt,dxs,dts); 


eeror the fluid 
Biot —(adts°2)/(kf*2*dxs*2); CC i — 2 a ATO), 
med = (2*dats*2)/(kf*°2*dxs*2) ; 


The matrices for the radiation boundary condition allowing the propagating 
modes to pass through the artificial boundary. 


oe dP 


anew = zeros(ct+l); 
for pm = O:st 


eure = Ebe(r+l,ct+1,pm,dxs,kf)-; 


anew = anew+acurr; 
end 
ml = (eye(ct+1) + (dts/(2*kf£*2))*anew); ml = inv(ml); 
m2 = (eye(ct+1l) - (dts/(2*kf*2)) *anew) ; 


me= trid(daxs,dts,ct+1,kf) ; 
% Step 3: vl and v3 calculated so as to be able to determine the amplitude 
vl Sel2(2*nt) ; 
v3 exp(i*mm*pi*x’) ; 


t Il 


+ Step 4: basel builds the I -beam pclor,pcoll ,coll,colr and fill 
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dP oP 


% 
% 


= 


a oo 


allow us to identify the boundaries of the solid domain and the 
corner nodes. 


{[b,rows, pcoll,pcolr,f1ill,coll,colr] = basen 2), 
Step 5: All values for the fluid and the solid are initially 


set to zero. 
un = zeros(size(m)); uc = un; uo = un; 


vn zeros(size(m)); vc vn; vo = vn; 


mn zeros(Ssize(m)); mc mn ;mold = mn; 


for k = linn 


Step 6: The elastic wave equation is applied to the interior 
of the solid 


un(2:n,2:n) = ces*uc(2:n,2.) )— UC(2-n, 2m)... 
+ rho(2) * (ve(lin=-17 2: nese isons hen) 
+ rho(1)*#(ue(2:n,-lin=)) + uc(2:- nest) ae 
+ rho(3)*(ve(lisn=-1,l:n=-1) + vets. 2s) 
- rho(3)* (ve Gen-1,32n+1) + vel3: nt) ae): 


Vn(2in,2:n) = Cest*tvcl(2: ny 23) = 902-2 se 

+ rho(1)*(veQ(len-T,Z2:ny & Vets: mee 2s eee 

+ rho(2)*(ve(2in,):n-1) + Vel2iny sn) 
+ rho(3)* (ucQGien=1, 12 n-1)) +e (32s nee. 
- rho(3)*(uc(1:n=-1,3:ntlje tees eee 


The periodic boundary conditions for the solid 


un(22n,1) = Cest*uc(2:n,;1) = ue (Zan) 2... 
+ rho(2)* (uc(itn=1,1) + velscnt) nn... 
+ rho(l)* (uc(2en, 2) + uc(2-n noe 
+ rho(3)*(ve(lin-1,n) + VeCent 2 ier 
= rho(3)*(ve(lin=1,2)"] ve(sent ran 


wn(22-hn neh) = un (2sn, 1) 
Vn(2in,1) = ces*vc(2:n, 1) —\voegzon 
+ rho(1)*#(ve(lin-1,1) +8 vc Gi) a 
+ Fho(2)* (ve(22n,2) + Ve( 2 nem)_) 
+ rho (3) * (uct ien-= 1, n e+ ie (Sie ee eee 
= rho(3)* (WeClin=1, 2) sues: nti. nee 
Vnt22n, n+l) = yen, 1) 


Step 7: The boundaries of the I beam and including the corner nodes 
are treated. 
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% 


fan ,vn} = bndf£nx (uc, vc, uo, vo, un, vn, rows, pcoll,cunx, Cvnx) ; 

{un ,vn} = bndfny(uc,vc,uo,vo,un,vn,rows,pcoll,pcolr,cuny,cvny) ; 
ain ,vn} = lcorny (uc, vc,uo, vo, un, vn,rows,pcoll,cuny,cvny) ; 

{un ,vn} = bndfpx(uc,vc,uo,vo,un,vn,rows,pcolr,cupx, Cvpx) ; 

{un ,vn} = bndfpy(uc,vc,uo,vo,un,vn,rows,pcoll,pcolr,cupy,cvpy) ; 
fine, vn} = leorpy (uc, vc,uo0,vo,un,vn,rows,pcoll,cupy,cvpy) ; 
{un,vn] EGOnrlo (uC, VC AG, VO, Un -vnenows ,pcoll, peolr rho) ; 


{un,vn} = t 


Step 8 : The 


Gomeq Ge, VeC,u0, VO,Un, Vn, rows, pcoll, pcolr,rho) ; 


solid and the fluid are coupled (Note the is where the forcing 


% function of the problem is contained ) 


all = me(l 
fun ,vn,mnp 


ole ble Me(2,. ) > 
djs— coup! (uc, vec,u0, vo, um, vapcupy,cvpy,k,dts,dxs,all,bll1,epss) ; 


% Spep 9: Both sides of the center spar for all values of u and v are zeroed 
% out so as to prevent pollution 

mimerows,cOll) = fill; vn(rows,coll) = £111; 

Weerows,COll) = fill; ve(rows,coll) = fill; 

Weerews,colr) = £111; vn(rows,colr) = fill; 

merows,coOlr) = fi11; vc(rows,colr) = fill; 

uo(rows,coll) = f111; vo(rows,coll) = fill; 

Worrows,cOlr) = fill; vo(rows,colr) = fill; 


ow 


Step 10: The scalar wave equation for the fluid 


* First the fluid/solid interface 
Mane: C) = rhof*(mnpd(2:c) + mc(2,2:c)... 
meee C- 1) + me(i,3:e+1)) +ccof*mce(1,2:c)... 
mola(),2:C) ; 


% and it’s periodic boundary condition 


eet) = rhnof*(mnpd(1,1) + mce(2,1)... 
fmeeiene)) + me(l,3)) +ccf*mc(1,1)..-.. 
emota(],1); 


macs, ctl) = mn(l,1); 


% The interior points of the fluid 
mueeemre 2. Cc) —]  rhol*(me(l:ir-1,2:¢) + me(3:rt+1,2:c)... 
mmeeeee ls C— 1) ot me(2:r,s:crl)) +tccf*mce(231r,2:c)... 
mrota(2:r,2:Cc) ; 


The radiation boundary condition 

Note: We are reguired to multiply a matrix by a row vector,but to do this 

it must be transposed to a column . Matlab takes the Hermitian transpose 

by default , so to ensure the correct signs we must void this effect by 

taking the conjugate of the transpose before we do our calculations. 

The process is then reversed so as the updated value has the correct dimension 


co oP of oO oP OW oO 
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$1 = con] (me) oe. 
$2 = con) (me(r+i 7. 
s3 = con} (moldiz+e..)):; 
int = m1l*(rcl*sl + t*s2 - m2*s3); 
mn(r+1,2) = conjunc 


% the periodic boundary condition for the artificial boundary 


mn(2:r,1) = rhof*(me(1sr—-l, 1) tene (2. tee 
mco(3:r+1, 1) +. me(2:5, Cc) )eteeerenc(2Z 1,1 es NO ae 


mn(2:r,Ctioy= Mab2:.,)) 

% Step 11: Calulates amd accumulates the values of the amplitude 
aa = abs((mn(r,:).*vl1)*v3); 

22.) — (22 aa | 


* Step 12: The values for u,v and m are updated. 


uo = uc; uc = un; 
vo = VC; ve = vn; 
mold = mc; mc = mn; 


end 
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function del2 - del funtion used for thr trapezoidal rule 

written by: Lt. Hugh Mc Bride USNR 

date ay -\ 0) a 2 

bulids a vector of appropiate length used to weight the elements of the 
quantity being integrated 


AP HO oP oP oO 


function d = del2(n) 

a = ones(1,n+1); 
more = .5; a(nt1) = .5; 
a2 = ones(1,nt+1)/n; 


d = a.*a2; 
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function lcorny - left corner facing in the negative y direction 
i.e unit normal = (0,-1) 

written by: Lt. Hugh Mc Bride USNR 

date ; APE 92 


lecorny performs the same function as bndfny except it only operates on 

4 points , those which lie at the extremities of the domain, the points 
as if fig( ) Periodicity is used for pseudonodes which lie outside the domain. 
This requires us to pick out each nodes individually (9 for each 4 points , 
making a total of 36 ) they are weighted in the same fashion 

as in bndfny and the updated values for u and v are calculated for both sides 


function (un ,vn) = leorny (uc, vc;,u0, vo,uUn, vn, Gows, peo eae, 


OP oP ao of CP oO OM CO OM oO 


oe of? 


un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values of v 

rows : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 

the nodes on the boundary. 

pcoll allows us to pick out the column location of those nodes 

to the left of centerline and using the periodicity we substitute 
this value into the corresponding point on the right of 
centerline. 


cuny; 
cvny; 


oe ol? 
ao 


Pout 2h) 
[12,32] 
[13,33] 


s$ize(pcoll); 
$1ze(rows); 
size(uc); 


the updated values of u to the left of centerline are 
calculated 


ucc = uc(rows(72)+1) peall (io 
ucecl = ue(rows (92) 4+2-peoll(i))). 
a 


ucr = uc(rows(j2)+1,pcoll(2 

ucl = uc(rows(j2)+1,33-1); 

ver = vce(rowsS(j32)4+1,peoll(2))- 

Vru = Ve(rows()2) +2, peel (257 

vlu = ve(rows(j}2)+2,33-1); 

vel = ve(rows(j2)+1,33-1); 

uol = uo(rows(}2)+1,pcoll(1)); 

ul = ¢c(1)*ucce = uol + €(2) *tecle+e(2)- (ices sue ere 


+0(4)*(vlu = vru) + €(S) *{ver— vel), 


ce) = Wells peo lina) 
UCUL ="NeE(2, peel. 
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Heri = ue(l,peoll(2)); 
Held = ye(1,)3-1); 


werd = vc(1,peoll(2)) 
wreout = ve(2,pcoll (2) 
vlul = ve(2,)33-1); 
Vell = vc(1,33-1); 


); 


Holl) = u0(1,pcoll(1)); 


eee tence! — uoll + c(2)*ucul +¢(3)*(ucrl + ucll)... 
+c(4)*(vlul - vrul) + c(5)*(verl- vell); 


% the updated values of v to the left of centerline are 
% calculated 


vuce = vc(rows(j2)+1,pcoll(1) 
vuccl = vce(rows(j2)+2,pcoll ( 
Wuer = vce(rows(j2)+1,pcoll ( 
mucl = vc(rows(j32)+1, 33-1 


~~ fA) ~~ 
et ee Ey 
ee 


VVer = uc(rows(j2)+1,pcoll(2)); 
eeeue= Uuc(rows(}2)+2,pcoll(2)); 
vvlu = uc(rows(j2)+2,33-1) ; 
vvel = uc(rows(j2)+1,j33-1); 


vuol = vo(rows(j2)+1,pcoll(1)); 


vul = d(1)*vucc - vuol + d(2)*vuccl + d(3)*(vucr + vucl)... 
+4(4)*(vvlu - vvru) + d(5)*(vvcr- vvcl); 


weer = vc(1,pcoll(1) 
mew = vc(2,pcoll(l 
wuerl = vo(1,pcoll ( 
vucll = vc(1,33-1 


); 
ee 


— A) — 


Meer! = uc(1,pcoll(2)); 

Meme = uc(2,pcoll(2)); 
vvlul = uc(2,33-1); 
vvecll = uc(1,33-1); 


ruoll = VO(1,peoll(1)); 

weil) *vucec] -— vuoll + da(2)*vucul +d(3)*(vucrl1 + vucll)... 
+d(4)*(vvlul - vvrul) + d(5)*(vvcerl- vvcll); 

% the values put in the correct positions 

% and since we have periodic boundary conditions, we insert the left 


% hand value into the correspondding right hand position 


Mme rOwsS(3}2)+1,pcoll(1)) = ul; 
Un(rows(j32)+1,33) = ul; 
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vn (rows(j2)+1,pcoll(1)) = vul; 


vn(rows(j2)+1,33) = vul; 
un(1,pcoll(1)) = ull; 
un(2,33)=suli: 


vn(1,pcoll(1)) = vull; 
Wn (1,33) = vue, 
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PO AO AP OO AP AO AP AP CP CP OP CW OP 


function lcorpy - left corner facing in the positive y direction 
i.e unit normal = (0,1) 

written by: Lt. Hugh Mc Bride USNR 

date >: Apr 92 


lcorpy performs the same function as bndfpy except it only operates on 

2 points , those which lie at the extremities of the domain, the points 
as in fig( ) Periodicity is used for pseudonodes which lie outside the domain. 
This requires us to pick out each nodes individually (9 for each 2 points , 
making a total of 18 ) they are weighted in the same fashion 

as in bndfpy and the updated values for u and v are calculated for both sides 


function {un ,vn} = lcorpy(uc,vc,uo,vo,un,vn,rows,pcoll,c,@q); 


PO AP AP AP AP OW DO oO OP CO 


of of # oO 


un : updated values of u vn :; updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values of v 

rows : used to identify the elements of the matrix which 

are zero for all times they also contain the row location of 

the nodes on the boundary. 

pcoll allows us to pick out the column location of those nodes 

to the left of centerline and using the periodicity we substitute 
this value into the corresponding point on the right of 
centerline. 


$c = cupy; 
% d= cCvpy; 


there are now only two values which need to be calculated 
as the other boundary facing in the positive 

y direction is at the fluid/solid interface and requires 
a special treatment 


{11,31} = size(pcoll); 
[i2,}2}] = size(rows) ; 
frogs) = $ize(uc); 


the updated u value 


Meen=— uc(rows(1)-1,pcoll(1)); 
ucu = uc(rows(1)-2,pcoll(1)); 
Weme— UC(rows(1)-1,pcoll(2)); 
ucl = uc(rows(1)-1, 33-1); 


mete VC (rows(1)-1,pcoll(2)); 
vru = vc(rows(1)-2,pcoll(2)); 
vlu = vc(rows(1)-2,33-1); 
vel = vce(rows(1)-1,33-1); 
Hol = uo(rows(1)-1,pcoll(1)); 


ul = c(1)*ucc - uol + c(2)*ucu +c(3)*(ucr + ucl)... 
+¢(4)*(vlu - vru) + c(5)*(vcr- vcl); 
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$ the updated v value 


vuce = ve(@ows (1) -1,;peet)(1)); 

vucu = vec(rows(1)-2,pcoll(1)); 
vucr = ve(rows(1)-1,pcoll(2)); 

vucl = ve(rows(1)-1,33-1); 


vver = uc(rows(1)—1, pcoll(2))) 
VVru = uc (rows (1) =2, peoll (2) 
vvlu = uc(rows(1)-2,33-1); 

vvel = uc(rows(1)-1, 33-1); 


i 


vuol = vo(rows(1)-1,pcoll(1)); 


vul = d(1)*vucec = vuol + d(2)*vucu +oaisi-(vuer +evucil). 
+0(4)*(vvlu - vvru) + d(5)*(vver- vvcl); 


$ using periodicity we update the values to the left and right of 
$ the center spar 


un(rows(1)-1,pcol1l(1)) = ul; 
un(rows(1)-1,33) = ul; 
vn (rows (1) <-1,pcoll (1) ) 
vn(rowsS(1)-1,j33) = vul; 


é 
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fumetion rbc - radiating bondary condition 
written by: Lt. Hugh Mc Bride USNR 
date : Mar 92 


rbc builds the matrix A(i,k) which is required by the radiation 
boundary condition 


AP APO AO AO oO OO 


fumecion Vv = rbc(j,1,pm,dxs,kf) 
% variables 
% j,1 -dimensions of the matrix 
%¢ pm - propagating modes 

a = zeros(j); 

fore! = 1:), 
por kK = i:l, 
a(i,k) = exp(sqrt(-1) *pm*pi*(i-k) *dxs) ; 


end 


end 


Gite, ) = .b%a(:,1) = a(:,K) = .5*a(:,K); 


v = sqrt(kf*2-(pm*pi) *2) *a; 


lel 


% function rho- del funtion used for thr trapezoidal rule 

% written by: Lt. Hugh Mc Bride USNR 

+ date >: Apr 92 

% generates a vector containing constants used repeatedly throughout 
% the program. 


function rho = rhov(kl,kt,dxs,dts) ; 


rhol =(dts%*2) /(k1*2*dxs%2) ; 


rnog 


(dts°2)/ (KtE°2*axs*2) ; 


rho3 


(dts*2/(4*dxs*2))*((1/k1*2)-(1/kt*2)); 


rho = [rhol rho2 rho3]; 


at 2 


% function stop - defines the stopping criteria for constructing the radiation 

% boundary condition i.e if the result of stop is 0 then only the fundamental mo 
% is allowed to propagate, 1 only the fundamental and the first mode are allowed 
% to prpogate and so on. 

% written by: Lt. Hugh Mc Bride USNR 

+ date : Mar 92 


function v = stop(k) 


% variables 

$k- 

$n - nth propagating mode 

n= 0; m= -1; 

while (k*2 - n*2*pi*2) > 0. 
m= m+1; n = nti; 


les 


OP OP OP OP OP OP OP OP GP OP OP dP oP 


OP dP OP dP AP oP OP GP AP CP CP 


Cf oo 


function tcorl13* -— treatment of corners ™“ianass 
written by: Lt. Hugh Mc Bride USNR 
date : “Apr 92 


tcorl13 applies the elastic wave equation at 

corners 1 and 3 as per fig( ) since they use the 

same difference formula. 

the corner nodes are located (1 first then 3) 

identified as p and q (pl ,qli in the case of corner 3) 

the necessary neighbours are picked off from the u and v 
Matrices and weighted accordingly and the new updated values 
for u and v are computed. 


function [un,vn] = tcorl13(uc,vc,u0,vo,un,Vvn, rows, pcoll, pcolr, rie 


variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : old values of u vo : old values of v 

rho : vector containing global constants 

rows : used to identify the elements of the matrix which 
are zero for all times they also contain the row location of 
the corner nodes. 

pclor and pcoll carries out the same function as rows for 
the columns, and the contain the column location of the 
corner nodes. 


(il, 92) 
[12,Jj2) 
[13,53] 


size(rows); 
size(pcolr); 
Size (pecoll), 


fo Wi 


we pick off the elments of rows and pcolr which 
identify the location of corner 1. 


rows(1il)<-1; 
peel (3) 


p 
q 


generate any required local constants 
ee = (2) — 2*rhotl])—2*rne(2)) 


the updated values of the corner nodes for u and v are 
calculated 


un(p,q) = cc*uc(p,q) - uo(p,q)... 
+ rho(1)*(uc(plq=)) 49uc(o, ae)... 
+ rho(2)* (uc(p=17q) + uc(otio a... 
~ 2*rho(3)*(vce(p-1,q-1) + ve(pti,q+l) + 2*vc(p,q))... 
+ 2*rho(3)*(vc(p,qtl) + ve(p,q-1) + vce(pt1,q) + vce(p-1,q)); 


vnip,q) = cc*vc(p,q) “= votp 7a). : 
+ Fhe 2) * (Ve( pigs) jae. Vee ncn ren 
+ rho(1)*(vce(p=-1,q) + ve(p+l,q))... 
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Zee eno(s) = (UeC(p=1,gq—-1) + uel(p+l, aii) 7.2*uc(p,q))... 
fee Eno (s) *(ue(p,qri) + we(p,q=-1)ms uc(prl,;d) + uc(p=1,q)); 


%$ the same process is repeated for corner 4 below 


rows (j1)+1; 
peor (12): 


pl 
foal 


Mimo, adil) =ce*uc(pi,gqi) - uo(pl,ql)... 
meant) =(UC(pl,gdli-1) + uc(pl,qitl))... 
meanoue *(UC(pl—1,ql) + uc(pl+1,ql))... 
mero ts) ~(VC(pi-l,qi-1) + vc(pl+l1,qit+1) + 2*vc(pl,gql))... 
+ 2%*rho(3)*(vc(pl,qiti) + vc(pl,qi-1) + ve(pi+1,ql) + vc(pl-1,ql)); 


wml, gl) =ce*vce(pl,gql) - vo(pl,ql)... 
Beene 2)*(vVC(pi,gi-1) + ve(pl,qlitl))... 
Papen.) *(vo(pl<1,q1) + vce(piti,ql)). 
- 2*rho(3)*(uc(pl- 1,q1-1)+uc(pl1+l1, qi+1)+2*uc(pl, 9/5) 
tee rio (5) = (UC( pl, Girl) rT uc(pl,gi-1) + uctpl+17qt) “ “uc(p1- 17 Gl); 
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DO AO AO AO AO AP GAO AO AO AO CM AO ao 


GO AO NO oP AP dP OP dP oO OP oP 


oe oe oe 


ae oe 


function tcor24 - treatment of corners 2 and 4 
written by: Lt. Hugh Mc Bride USNR 
date : Apr 92 


tcor24 applies the elastic wave equation at 

corners 2 and 4 as per fig( ) since they use the 

same difference formula. 

the corner nodes are located (2 first then 4) 

identified as p and gq (pl ,ql in the case of corner 4) 

the necessary neighbours are picked off from the u and v 
matrices and weighted accordingly and the new updated values 
for u and v are computed. 


function (un,vn)} = tcor24(uc,vc,uo,vo,un,vn,rows,pcoll,pcolr,rnae 


variables 

un : updated values of u vn : updated values of v 

uc : current values of u ve : current values of v 

uo : Old values of u vo : old values of v 

rho : vector containing global constants 

rows : used to identify the elements of the matrix which 
are zero for all times they also contain the row location of 
the corner nodes. 

pclor and pcoll carries out the same function as rows for 
the columns, and the contain the column location of the 
corner nodes. 


size(rows); 
size(pcolr); 
size(pcoll); 


Celah ais 
(12,J2] 
[13,J3) 


ft od 


we pick off the elments of rows and pcolr which 
identify the location of corner 2. 


p = rows(11)-1; 
q = pcolrii2) 


generate any required local constants 


6e = (2° = 2%Eho(1) -2*uho(2jye 


the updated values of the corner nodes for u and v are 
calculated 


un (p,q) = cc*uc(p,q) = US(p, a). 
+ rho()) *(Uuci(p, 4-1) ue (eGo aa 
+ rho( 2) * (Uc(p=1 7d) -avuetp. ld) )-..- 
+ 2*rho(3)*(ve(p+1,q-1)) +5 ve(p-1,q+1) + 2*ve(p,q)).-- 
- 2*rho(3)*(vc(p,qt1l) + vc(p,q-1) + vc(p+l,q) + ve(p-1,q)); 


vn(p,q) = cc*vc(p,q) - vo(p,q).-.-. 
+ rho(2) *(vel(p,q-1) + VeGera tiie... 
+ Eno(1) *(ve(paig) + VveGgr- 1 q)) =... 
+ 2*rho(3)*(uc(prl, g-]je uci p-1,q+1)) +e2%Ue (p,q). 
- 2*rho(3) *(uc(p,q+1) + uc(p,q-1) + uc(p+l1,q) + uc(p=1,4q)); 


1d6 


%$ the same process is repeated for corner 4 below 


pl = rows(j1)+1; 
ql= pcoll(3j3); 


Bma(pl,ql) = cce*uc(pl,gql) - uo(pl,ql)... 
penne l)*(uc(pl,qi-1) + uc(pl,qli+1))... 
feo (2) *(uc(pl-1,ql) + uc(pl+l1,qil)).. 
+ 2*rho(3)*(vce(pl+l, qi-1) + vce(pl-1, qi+1) te 22 Ve Cp, G1 ))). 
=e2ernota)~(ve(pl,qitl) + ve(pl,gqi-1) + ve(pl+1,qli) + ve(p1-1,q1)); 


Vimopl,ql) = ce*ve(pl,ql) - vo(pl,ql)... 
Teene,2)*(VC(pl,gqi-1) + vce(pl1l,qitl1))... 
pernotl)*(ve(pl-1,q1) + vec(pl+l1,ql))... 
peemrne a)*(UC(pl+1,gqi-1) + uc(pl=1,q1l+1) + 2*uc(pl,ql))... 
ened )(UC( ol git ll) teuc( pl, gi=1) + uc(plt+1,ql) + uc(pi-1,q1)); 
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function trid - tridiagonal 
written by: Lt. Hugh Mc Bride USNR 
date : Mar 92 


trid generates a tridiagonal matrix for the 

radiation boundary condition of the fluid 

the elements of the sub and super diagonal are (dts/kf*dxs) *2 
the main diagonal component is 2 - 4(dts/kf*dxs) *2 

the (n,2) and (1,n-1) contain (dts/kf*dxs)*2 to satisfy the 
periodic boundary conditions. 


JP dP AP AP AP GP dP AP de dP 


function m = trid(dxs,dts,n,kf) 


% variables 

$ dxs : scaled spacing 

$ dts : scaled time step 
$n : dimension of matrix 
* kf : scaled constant 


$ an identity matrix for the elements of the main diagonal 


dl = eye(n); 


oo 


the sub and super diagonals 
dad2 = diag(ones(n-1,1),1) + diag(ones(n-1,1),-1); 
$ the required co-efficients 
rho = dts/(kf*dxs); rho = (rho*2); 
%$ generates the required matrix 
dad =2*dl -4*rho*dl +rho*d2; 


d(1,n-1)= rho ; d(n,2) = rho; 
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function ucnx - u coefficients/constants for the boundary facing the 
negative x direction. 

written by: Lt. Hugh Mc Bride USNR 

date : Apr 92 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(-1,0) for the corresponding u 
and v values. 


JOC AC AC HO AO AO dO GO 


function cunx = ucnx(kl1,kt,dxs,dts) 
% variables 
% kl - scaled longitudinal speed 
* kt - scaled transverse speed 
% dxs - scaled spacing 
%* dts - scaled time step 
% the coefficients are calculated and stored in the vector cunx 
% for use in bndfnx 
Cl = 2 - 2*(dts*2/dxs*2) *((1/k1*2)+(1/kt*2)); 
e2)— 2*dts*2/ (adxs*2*kl*2); 
c3 = (dts*2)/(dadxs*2*kt%*2); 
e4 = ~(dts*2/(2*dxs*2) )*((1/k1*2)-(1/Kt%*2)); 
c5 = (dts*2/(2*dxs*2))*((1/k1*2)-(3/kt%*2)); 


cunx = [cl c2 c3 c4 c5]j; 
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function ucny - u coefficients/constants for the boundary facing the 
negative y direction. 

written by: Lt. Hugh Mc Bride USNR 

date : Apr 92 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(0,-1) for the corresponding u 
and v values. 


AP AP AP AP AP oP oO AO 


function cuny = ucny(kl,kt,dxs,dts) 


variables 

kl - scaled longitudinal speed 
kt - scaled transverse speed 
adxs - scaled spacing 

dts - scaled time step 


OP AP AP ah ah 


the coefficients are calculated and stored in the vector cuny 
for use in bndfny 


dP oo 


Cl = 2°] "2% (dts72/dxs*2) * (1 /ki@ zee ke 2); 


C2 = 2*(dts*2/(dxs*2*kt%*2) ); 


9) 
Ww 
Il 


(dts*2) /(dxs*2*k1%*2) ; 
c4 = -(dts*2/ (2*dxs*2) )*((1/k1%2)-(1/kt*2)); 
c5 = (dts*2/(2*dxs%*2) )*((3/kt*2)-(1/k1*2)); 


cuny = [cl ¢c2 C3 7e4y coir, 
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$ function ucpx - u coefficients/constants for the boundary facing the 

$ positive x direction. 

% written by: Lt. Hugh Mc Bride USNR 

% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal(1,0) for the corresponding u 
%$ and v values. 

% 
f 


unction cupx = ucpx(kl,kt,dxs,dts) 
% variables 

%* kl - scaled longitudinal speed 

+ kt - scaled transverse speed 

% dxs - scaled spacing 

% dts - scaled time step 

% 
$ 


the coefficients are calculated and stored in the vector cupx 
for use in bndfpx 
Cl = 2 - 2*(dts*2/dxs*2) *((1/k1*2)+(1/kt*2)); 

ee. — 2*(dts°2/ (dxs*°2*k1l*2) ); 


faitss2) / (Axss2*kt~2) ; 


Q 
Ww 
II 


Bee— (dts-2/ (2*dxs~2) ) *((1/k1*2)-(1/kt*2)); 
wee (ACS 2/ (2 *axXS>2)) = ( (1 / Kl 2) —(3/kKeE°2)); 


Supe — (Cl c2 c3 c4 ¢5]j; 


al 


function ucny - u coefficients/constants for the boundary facing the 
positive y direction. 

written by: Lt. Hugh Mc Bride USNR 

date $ Mpr 92 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(0,1) for the corresponding u 
and v values. 


OP dO AO dO HO AO oO CO 


function cupy = ucpy(kl,kt,dxs,dts) 


variables 

kl - scaled longitudinal speed 
kt - scaled transverse speed 
adxs - scaled spacing 

dts - scaled time step 


JP fm dO oO dO 


% the coefficients are calculated and stored in the vector cupy 
% for use in bndfpy 

cl = 2 = 2%* (ats 27axs" 2) *(€1/K) 2) ree 

c2 = 2*(dts*2/(dxs*2*kt*2) ); 


c3 = (dts*2) /(dxs*2*kl*2); 
c4 = (dts*2/(2*dxs*2))*((1/k1*2)-(1/kt*2)); 
c5 = (dts*2/(2*dxs*2) )*((1/k1*2)-(3/kt*2)); 


cupy = (cl ¢€2°%e3 C4 eo); 
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function venx - v coefficients/constants for the boundary facing the 
negative x direction. 

written by: Lt. Hugh Mc Bride USNR 

date : Apr 92 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(-1,0) for the corresponding u 
and v values. 


OP AP AP AO HO AO dO HO 


function cvnx = vcnx(kl,kt,dxs,dts) 


variables 

kl - scaled longitudinal speed 
kt - scaled transverse speed 
adxs - scaled spacing 

dts - scaled time step 


AO dP AO AO HO 


the coefficients are calculated and stored in the vector cvnx 
for use in bndfnx 


oe 


di = 2 - 2%(dts*2/dxs*2) *((1/k1*2)+(1/kt*2)); 
d2 = 2*(dts*2/(dxs*2*kt*2)); 
a3 = (dts*2) /(dxs*2*kl*2); 
a4 = -(dts*2/(2*dxs*2)) *((1/k1*2) -(1/kt*2)); 
a5 = (dts*2/(2*dxs*2))*((3/kt*2)-(1/k1*2)); 


cvnx = [dl d2 d3 d4 d5]; 
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function vcny - v coefficients/constants for the boundary facing the 
negative y direction. 

written by: Lt. Hugh Mc Bride USNR 

date ; Apr 932 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(0,-1) for the corresponding u 
and v values. 


AO AO HO AO AO AO oP oP 


function cvny = vceny(kl,kt,dxs,dts) 
% variables 
% kl - scaled longitudinal speed 
%* kt - scaled transverse speed 
% dxs - scaled spacing 
% dts - scaled time step 
% the coefficients are calculated and stored in the vector cvny 
% for use in bndfny 
dl = 2 = 2*(dts*2/dxs*2) *((1/k1*2)+(1/kt*2) ); 

d2 = 2*(dts*2/(dxs*2*k1%2)); 

a3 = (dts*2) / (dxs*2*ke72); 
qd4 = -(dts*2/ (2*dxs*2) )*((1/k1%2) -(1/kt%2) ); 
dS = (dts*2/(2*dxs*2))) (4 kt oye tae 


cvny = [dl d2 d3 a4 a5); 
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% function vcpx - v coefficients/constants for the boundary facing the 

% positive x direction. 

%$ written by: Lt. Hugh Mc Bride USNR 

% date : Apr 92 

% provides the coefficients developed by applying the Ilan and Lowenthal 
% technique to the boundary with unit normal(1,0) for the corresponding u 
% and v values. 

% 
£ 


unction cvpx = vcpx(kl,kt,dxs,dts) 


% variables 

% kl - scaled longitudinal speed 
% kt - scaled transverse speed 

% dxs - scaled spacing 

% dts - scaled time step 


the coefficients are calculated and stored in the vector cvpx 
for use in bndfpx 


JP oO 


dl = 2 -2*(dts*2/dxs*2) *((1/k1*2)+(1/kt*2));3 
d2 = 2*(dts*2/(dxs*2*kt*2)); 
a3 = (dts*2)/(axs*2*k1*2) ; 
pe dts*2/(2edxs-2) )*((aykl- 2) —(1/kt-2) 
d5 = (dts*2/(2*dxs*2))*((1/k1*2)-(3/kt*2)); 


Cvpx = [dl a2 d3 d4 45}; 


15 


function vcny - v coefficients/constants for the boundary facing the 
negative y direction. 

written by: Lt. Hugh Mc Bride USNR 

date Apr SZ 

provides the coefficients developed by applying the Ilan and Lowenthal 
technique to the boundary with unit normal(0,1) for the corresponding u 
and v values. 


AP AP A HO AP A OP DO 


function cvpy = vepy(kl,kt,dxs,dts) 
variables 

kl - scaled longitudinal speed 
kt - scaled transverse speed 

adxs - scaled spacing 

dts - scaled time step 


A AP A oO oO 


the coefficients are calculated and stored in the vector cvpy 
for use in bndfpy 


ow 


Gl = 2 - 2*(dts*2/dxs"2)*((1/kl(2) 4 Ge hee 


QO, 
N 
Il 


2* (dts*2/(daxs*2*k1*2)); 


QO, 
W 
ll 


(dts*2)/(axs*24*ktu- 2) 
a4 = (dts*2/(2*dxs*2))*( (7k) 2) Oy ke 2) 
d5 = -(dts*2/ (2%*axs°2)) * ( (17 Kile 2) = xe 2, 


cvpy = {dl d2 d3 a4 a5, 
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